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ABSTRACT 


The objective of this research is to develop the theory and 
associated numerical technique for the iterative design improvement 
of the compensation for linear, time- invariant control systems with 
multiple inputs and multiple outputs. The multivariable capabili- 
ties allow system suboptimization of several control loops with 
coupled characteristics. A strict constraint algorithm is used 
in obtaining a solution of the specified constraints of the control 
design. The result of the research effort is the Multiple Input, 
Multiple Output Compensator Improvement Program (CIP) . 

The objective of the Compensator Improvement Program is to 
modify in an iterative manner the free parameters of the dynamic 
compensation matrix so that the system satisfies frequency domain 
specifications. In this exposition, the underlying principles of 
the multivariable CIP algorithm are presented and the practical 
utility of the program is illustrated with space vehicle related 
examples. Further, the capabilities of and possible extensions to 
the algorithm are delineated. 
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CHAPTER I 


INTRODUCTION AND SURVEY OF PREVIOUS WORK 
Introduction 

The space age era has challenged the control theorist to 
examine the fastest and most efficient means of system design and 
analysis. In this regard, increasing emphasis has been stressed on 
the utilization of digital computers in modern control theory. Much 
has been written on employing digital computer control methods for 
single input systems, but the complexity of systems may require more 
than this simple approach — that is, for a better control the sensing 
of many system parameters and inputs as well as their relationships 
to one another must determine the control law. The facet of control 
theory requires exploiting multiple input, multiple output tech- 
niques. Further, the problem of producing the ’best' compensated 
system becomes one of satisfying physical restrictions of system 
parameters as well as digital computer limitations. 

Modem trends in engineering systems are toward greater com- 
plexity, due mainly to the requirements of complex tasks and the 
necessity for accuracy. Sophisticated technology involves systems 
that are described adequately only by numerous variables; thus, 
these high-order complex systems are responsible for the dichotomy 
of the academic and industrial attitudes. The naivete'' of low-order 
system textbook orientation must be abondoned, and with the essence 



2 


of the classical methods espoused, analysis of large-scale systems 
made possible with modem methods. The necessity of meeting 
increasingly stringent requirements on the performance of control 
systems, the increase in system complexity, and the accessibility 
of large-scale computers has forced modem theorists to reexamine 
the problem of interaction of multiple control inputs. A complex 
system may have many inputs and outputs which may be interrelated in 
a complicated manner. To analyze such a system, it is essential to 
reduce the complexity of the mathematical expressions as well as to 
resort to computer algorithms for the tedious computations. 

Many process control systems have multiple inputs and outputs. 
Generally, there is no assurance that changes in one reference input 
will affect only one output; thus, the inputs and outputs are not 
decoupled but interact with one another. Analysis of such inter- 
actions of all inputs with all outputs is a difficult task, particu- 
larly by classical means. It is for this reason that the investi- 
gation and development of a computerized algorithm for compensation 
of multiple input, multiple output systems with interactions of 
parameters are so important. 

In conventional control theory, generally only the input, out- 
put, and error signals are of concern; the design and analysis are 
developed using transfer functions, together with a variety of 
classical techniques such as root locus, Nyquist, Bode, etc. The 
appealing characteristic of conventional control theory is that it 
is based on the input-output relationships of the system; that is. 


the transfer function. 
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The main disadvantage of conventional theory is that, primarily, 
it is applicable only to linear time-invariant systems having a 
single input and output. It is difficult to apply to time -varying 
systems, nonlinear systems except the simplest cases, and to multi- 
ple input, multiple output systems. Thus classical techniques are 
not amenable generally to the design of optimal or adaptive con- 
trollers . 

With the advent and abundant utilization of digital computers, 
computer-aided control system design has become a popular research 
topic. However, most contemporists abandoned the classical tech- 
niques and explored the realm of optimal control theory delving in 
systems described by state variables, i.e., a set of first order 
differential equations with design- objectives described in a cost 
functional. Although many enlightening methods and results have 
been evidenced, the weaknesses are inherent. For example, the 
optimal control law iS' extremely dependent upon the proposed cost 
function and in many cases the correct cost function debatable [1]? 
Furthermore, all states are assimied available for feedback. Even 
with observer theory to reproduce iinmeasured variables, and subse- 
quently allow fewer measurements, the computer storage required is 
excessive. 

The current impression, that in order to utilize computer 

facilities the control problem must be implemented .in state variable 

* * , 

form, must be resolved'. There is, generally speaking, no substitute 
for state variable techniques when applied to simple control 


*Numbers in square brackets designate referenced items. 
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systems, but when solving relatively large plants the handling and 
storing of state matrices can be unthinkable. Thus if the classical 
theory can be extended to include modern concepts , as well as com- 
puter techniques, then possibly a more effective control tool will 
evolve. An interesting aspect of this approach is that it clearly 
shows the return to prominence of the classical frequency domain 
techniques in modem system analysis. 

Literature Survey of Previous Work 

In the past, several papers have appeared discussing various 
approaches to computerized classical design of control systems. 
Generally, performance specifications are satisfied by frequency 
response and root locus methods using trial-and-error procedures. 
According to [2], there are three digital techniques which appear to 
have merit in this regard: Automatic Frequency Domain Synthesis of 

Multiloop Control Systems (AUTO) , Compensator Improvement Program 
(CIP) , and Computerized Optimization of Elastic Booster Autopilots 
(COEBRA) In addition to these, .developments by Nail on the Eigen- 
value Encouragement Technique , [3] , as well as. .Mancini’s Computer 
Aided Control System Design Using Frequency Domain Specifications 
(CALICO) [4] and Vines’ Computer Automated Design of Systems (CADS) 
[5] are also of interest. Each of these algorithms is summarized, 
and the disadvantages and limitations are discussed. 
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Automatic Frequency Domain Synthesis of 
Multiloop Control Systems (AUTO) 

The algorithm AUTO [6], developed at the Aerospace Corporations 
was designed to aid in the synthesis of compensation for multi- 
loop, time-invariant control systems exemplified in Figure 1. In 
this figure, the plant P(s) is assumed to have fixed characteris- 
tics, a' single control input, and multiple outputs. With the feed- 
back path broken at a an open-loop transfer function is defined as; 

c(s) = C(s)/R(s) (1.1) 

assuming R(s) is unity for all frequencies s = j( 0 , then C(jto) is 
the open-loop frequency response. 

The philosophy of AUTO is to fit the open-loop frequency 
response C(jw) to the desired open-loop response C(j(jo) by selection 
of the parameters of the compensators, G (s), G (s) , ..., G (s) . 

X ^ X"! 

The compensators are varied algebraically by making incremental 
changes in their parameters. Each compensator is assumed to be a 
rational function of the form: 

M . N 

^ I a^^(s)^“V[l + I 3 (1.2) 

i=l i=2 

where M and N are chosen by the designer and a , b represent 

iCX K>x_ 

the compensator coefficients. The number of compensators is 
dictated by the number of plant outputs. 



PLANT 

DYNAMICS MEASURED STATES COMPENSATION 



R(s) IS THE SCALAR INPUT RESPONSE 
C( s ) IS THE OUTPUT RESPONSE 
P(s) IS THE H X 1 MATRIX TRANSFER 
FUNCTION DESCRIBING THE PLANT 
Xj(s) REPRESENTS THE MEASURED STATE 
OF THE (J)TH CHANNEL; 

J = 1, 2, ... , H 

Gj(s) REPRESENTS THE (J)TH CHANNEL 
COMPENSATION; J = 1, 2, ...,M 


Figure 1. General Configuration for a Single 
Control Input System. 
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The measurement of closeness J between the actual response C(s) 

y\ 

and desired response C(s) is the mean square difference at a set of 
selected frequency points k = 1,2,...., K, that is 

J = I 1(6* - C*)^ W^W(C - C) 1 I (1.3) 

where the asterisk (*) denotes complex conjugate and 

c'^ = [C(jO)^) C(jW2) ... C(jWj,)] . . (1.4) 

W is a diagonal matrix of the form 

W 0 

0 ¥ 

0 0 

that is -used t o as sign different weights to the errors of different 
frequency points. The compensators are designed hy varying their 
coefficients so that J is minimized by a gradient search method. 

The directional vector along which the search is made is the gradient 
of J with respect to percentage changes in compensator coefficients; 
this type of directional vector is used to avoid premature con- 
vergence from encountering steep valleys or ridges associated with 
the function J. 

Compensator Improvement Program (CIP) 

CIP is a computerized design algorithm for aiding in the com- 
pensation synthesis for multiloop, time-invariant control systems 
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of the form of Figure 1. The basis of this algorithm is that with 
the feedback loop broken at a, the compensation design is accom- 
plished by satisfying certain frequency response specifications on 
the open-loop response C(jw)/R( 3 to) . CIP design specifications 
include the capability of obtaining gain margins (GM) , phase mar- 
gins (PM) , stability margins (SM) , and attenuation margins (AM) . 

Both the gain and phase margins have the normal definitions except 
their measurements in CIP are converted to distances from the 
(-1 + jO) point in the GH(ju) plane. The stability and attenuation 
margins are defined as follows: [7] 

Definition I 

For a closed— loop stable system whose open-loop fre— 
j quency response is described ,by GH(jw)., a stability 
' margin (SM) is defined as a relative minima of the real 
function^ 

Jl.O +, GH(jO)) j 

Definition 2 

An attenuation margin (AM) of the GH(j(o) frequency 
response for a band of frequencies such that lOj < w < Oi^ 
is defined as a relative maxima of the real function, 

|GH(jw)P, when we (w^jW^). 

Gain, phase, and, stability margins establish desirable amounts 
of phase stabilization; whereas, the attenuation margin is used to 
insure proper amounts of gain stabilization. CIP was developed with 
the objective of improving the frequency response from iteration to 
iteration. Two possible modes of operation are available: the Sum 

Improved Frequency Response (SIFR) , and the Total Improved Frequency 
Response (TIFR) . The user must select one of these modes to control 
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the amount of the incremental changes made in the compensator 
parameters in initiating the program.. Generally the SIFR mode 
allows coarser changes than the TIFR mode. 

GIF employs the mathematical programming tool the Constraint 
Improvement Algorithm (CIA) . This algorithm possesses the unique 
capability of producing a directional change vector for the compen- 
sator coefficients that insures the existence of a Total and/or Sum 
Improved Frequency Response. 

Computerized Optimization of 
Elastic Booster Autopilots (COEBRA) 

COEBRA design is achieved by solving a sequence of constrained 
optimization problems by minimizing a cost function. The cost 
function, in terms of frequency response of time domain specifica- 
tion's, is subject to a set of inequality constraints. The frequency 
response specifications include the classical phase and gain margins 
whereas, the angle of attack is included in the time domain specifi- 
cations . 

COEBRA employs the linear programming tool, the Simplex 
Algorithm, to obtain a solution; whereas, the design problem for 
which COEBRA was developed is. nonlinear in nature. However, if the 
cost and constraint function are approximated by a truncated Taylor 
series expansion, the problem becomes a linear one and a solution 
is obtained through a parametric programming procedure. . In ob- 
taining the truncated Taylor series expansion a finite difference 
technique yields the necessary partial derivatives. 
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Eigenvalue Encouragement Technique (POP) 

POP, a numerical technique for the iterative design of linear, 
time-invariant control systems, attempts to design dynamic feedback 
compensation by affecting the closed-loop eigenvalues in a desirable 
manner. This technique encourages the eigenvalues to migrate either 
toward or further into the left half plane, or toward other speci- 
fied values. This encouragement process is accomplished by solving 
an unconstrained minimization problem with a selected cost function. 

The algorithm is based on Danilevskii' s method of generating a 
characteristic polynomial and the assumption that the compensation is 
dynamic feedback. A unique relationship between a determinant and 
the partial derivative operation is applied to the system character- 
istic |Xl - A I ; the result is the partial of the closed- loop 
eigenvalues with respect to parameters in terms of 2p(n+l) determi- 
nants where n is the order of the plant and compensation A matrix, 
p the compensator order, and X the eigenvalue. By determinant 
manipulations the necessary 2p determinants are evaluated by 
Danilevskii* s methods yielding results for all the eigenvalues. 

Computer Aided Control System Design Using 
Frequency Domain Specifications 

CALICO, the computer-aided compensator design algorithm by 
Mancini, utilizes the constrained optimization method introduced by 
M. J. Box. This technique requires the desired open- loop frequency 
response be specified for discrete frequency points. The minimi- 
zation routine varies the compensator parameters in such a manner 
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as to minimize a cost functional based on the difference between the 
actual and desired frequency response of the compensated system* 

With the desired response as input, the author incorporates the 
normal frequency domain specifications such as gain and phase margins 
into the overall cost functional each time the algorithm evaluates 
the frequency response of the open-loop system, and thus eliminates 
the need for specialized computations to determine margin satis- 
faction. 


Computer Automated Design of Systems 

The automated digital computer technique by Vines, CADS, is a 
control system compensator design oriented in the time domain. In 
order to minimize a specific cost functional and set the free system 
parameters, the technique requires as input the desired output 
response and system description. The minimization technique BOXPLX 
by M. J, Box is employed. To simulate the system to be optimized, 
the author chose commonly used transfer functions which were reduced 
to first order linear differential equations. The equations are 
programmed so that the transfer function blocks can be cascaded by 
data card input. Several nonlinear transfer blocks are also availa- 
ble. The program simulates the system with known parameters and 
then allows all free parameters to be fixed by the optimization 
routine in achieving the desired response. 
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A Comparison 

Iii summary, each of these algorithms has advantages and 
limitations. With the exception of COEBRA and POP, these methods 
completely ignore the possibility of multiple inputs and their 
interactions with system parameters. Unfortunately, POP and COEBRA 
have several theoretical and computational limitations. 

For example, POP is a numerical technique that attempts to 
minimize a cost function composed of 'soft' constraints. This method 
should minimize the cost function, but in a practical sense, (in 
terms of relative stability, etc.) the final system may not be any 
better than the original. Further, the practical limitation of 
computer storage and run time may prove the infeasibility of apply- 
ing this method to large systems. In fact the run time is approxi- 
mately proportional to n^ where n is the system order including 
compensation, and systems above 40th order require more than 128K 
words of core storage on a UNIVAC-1108. Another unfortunate 
obstacle of this technique is the inherent problem of relating fre- 
quency domain design specifications such as phase and gain margins 
to closed-loop pole locations of large systems. 

COEBRA minimizes a cost function subject to a set of inequality 
constraints; the cost function is used to optimize gain and phase 
margins, rise time, percent overshoot, etc. , while the constraints 
insure that the performance measurements do not degrade from 
iteration to iteration. The directional vector is determined so as 
to minimize the cost function and not violate the constraints. 
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Thus, COEBRA also possesses the property that the final design will 
not be worse than the initial. COEBRA employs a method of finite 
differences^ to determine partial derivatives. This introduces 
numerical inaccuracies that jeopardize the practical utility of this 
program. From a user's view, COEBRA is difficult to enable; it 
requires a thorough knowledge of the programming techniques for a 
user to achieve a useful design. COEBRA also requires excessive 
computer core storage and run time; on the UNIVAC-llOS, 66K words of 
storage must be available; whereas, both AUTO and CIP can be 
executed in less than 32K words for systems of equivalent order [2]. 
Perhaps it was in this regard that the author found it necessary to 
restrict the design to no more than eight bending and/or slosh modes. 

AUTO assumes a single control input plant described by fre- 
quency response information in the form of complex numbers. Although 
AUTO’ appears easy to use, the judicious selection of weighting con- 
stants for every frequency component is a designer's nightmare; in 
some instances, no choice of constants will yield the desired design 
specifications. Like POP, AUTO seeks to minimize a cost function 
composed of soft constraints by a gradient optimization technique; 
the weighted mean square difference between the actual and desired 
frequency response indicates the measure of closeness to the desired 
results. 

CADS, as a time domain method, accepts only first order dif- 
ferential equations in describing the compensation. The computer 
time and storage are a function of the system order and the search 
area on the upper and lower bounds of the system parameters. The 
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routine requires a good guess on the original parameters in order 
to obtain a workable solution. According to Vines, the optimiza- 
tion routine BOXPLX may continue to search for a reduced cost 
functional to within- some significant digit even though a practical 
solution already has been determined. 

CALICO, although frequency domain oriented, is similar to CADS 
in its applicability to single input-output systems and its use of 
the BOXPLX optimization routine. Again the technique is cursed with 
an optimization method that becomes more inefficient and time con- 
suming as the system parameters increase, 

CIP has the limitation of applicability to a single input, 
multiple output system. In addition,- CIP requires much data in the 
form of frequency response information. Unlike the aforementioned 
method, however, CIP is not an optimization technique and does not 
attempt to maximize or minimize a cost function; rather, CIP searches 
for a ’ subop timal* feasible solution by satisfying a set of strict 
constraints that measure the performance of a design. The computer 
run time is proportional to the number of frequency points used to 
describe the plant. If properly used, CIP results in a final design 
that will be better than the initial. 

CIP has proven its merit in the service of space vehicle con- 
trol according to NASA contract reports [7,8]. If CIP could be 
extended to handle designs for plants with multiple inputs and 
multiple outputs, it obviously would be superior to any of the 
aforementioned methods. Further, the Inclusion of a multiple input. 



15 


multiple output capability would improve greatly the utility of this 
technique as a design aid. 

The Research Area 

The objective of this research is to achieve compensation for 
a multiple input, multiple output control system by developing an 
algorithm to facilitate fast and practical compensator design with 
maximum computer economy while minimizing designer effort. The 
relative stability method of the Constraint Improvement Algorithm 
[7] by McDaniel and Mitchell has been chosen to detemine system 
perfoirmance specifications by frequency domain techniques. The work 
presented in this exposition develops the theory and associated 
modifications necessary to extend the CIP type algorithm to the 
multivariable control system. With permission of the author of the 
original CIP [7] this program will now be known mnemo technic ally 
as CIP, since much of the philosophy and many of the techniques of 
the original algorithm have been retained and extended. 
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DESCRIPTION OF THE DESIGN ALGORITHM 
An Overview of the Design Problem 

Figure 2 is a schematic representation of a multivariable, 
linear, time-invariant feedback control system. This multivariable 
system may be viewed as n coupled feedback systems — one 'for each 
element of the input vector. The loop transfer function for the 
kth system is obtained by opening the feedback 'path at and then 
determining the response C^(s)/R^(s) with all other input R's set to 
zero . 

With this view of the multivariable feedback system, the de- 
signer is faced with the problem of synthesizing controllers of n 
interacting systems. ^TJsing classical feedback theory a controller 
may be designed so that the open-loop frequency response satisfies 
a set of frequency response design objectives. In theory this 
approach easily is extended to the multivariable system. However, 
in this case simultaneous designs of the controllers must be made so 
that the n open-loop frequency responses satisfy n. sets of design 
objectives. The simultaneity of the designs is required because of 
the implicit functional relationship between the design objectives 
of the individual systems; e.g,, a controller may affect the open- 
loop frequency response of one system in a desirable manner, while 
adversely affecting the response of another system. 
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Using classical frequency response techniques, the design of 
a control system to satisfy a few objectives can be acconq>lished 
with manual calculations. However, as the complexity of the system 
and the nximber of design objectives increase the development of the 
controller requires the aid of a high-speed digital computer. 

In order to use efficiently the digital computer in a design 
capacity it is necessary to have a design algorithm that is amenable 
to digital computation; in general, such algorithms are iterative in 
nature. The Compensation Improvement Program [7] is an algorithm of 
this type that has been developed to facilitate in the design of 
controllers for the class of systems of Figure 1. With the loop 



j is the controller index (j = 1,2,..^,M),, so that the open-loop 
frequency response, C(s)/R(s), satisfies specified requirements. In 
this study a design algorithm is developed to facilitate the design 
of controllers for multivariable feedback systems. 

The Algorithm Design Philosophy 

In order to accomplish logically the design algorithm the compu- 
tational flow diagram of Figure 3 Has been developed. The descrip- 
tion of the multivariable configuration requires discrete frequency 
data from each input to each output; whereas, the initial compen- 
sation for each controller input is described by a matrix of trans- 
fer functions. With this information the open-loop frequency 
response is obtained for each of the n coupled systems by deter- 
mining the associated subsystem response with one loop open at a 
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Figure 3. Computational Flow Diagram of the CIP Algorithm. 
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time. Likewise, for each subsystem, a set of critical points, that 
isj frequencies at which margins of stability or attenuation occur, 
is determined. It is possible to demand any number of margin 
requirements for gain, phase, stability, and attenuation radii. 
Further, these margins can be manipulated so as to make the con- 
straints or specifications frequency dependent. An active list of 
radii requirements is then prepared to alleviate any margins al- 
ready satisfied. From the open-loop frequency response data and the 
compensator coefficients, the gradient vectors of the active margins 
for each subsystem are calculated. Using these gradient vectors, a 
directional vector that can yield improvements in all active margins 
is determined. The free compensator coefficients are varied then in 
accordance with this directional vector. From this design the total 
response is checked to determine any margin radii not satisfied, arid 
the process continues in an iterative manner tmtil all specifications 
are met or user control, such as maximum computer time or iterations, 
forces a stop. 

Theoretical Concepts Associated with the Algorithm 

Calculation of Compensated Open-loop Frequency Response 

In order to develop a CIP type algorithm for designing the 
controller for the multivariable system, it is necessary to have 
equations for calculating the open-loop frequency response for each 
subsystem and equations for calculating the change in each objec- 
tive function (performance measurements) with respect to variations 
in the free parameters of the compensation. First attention is 
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focused on the calculation of open-loop frequency response infor- 
mation. From Figure 4 the output frequency response [C(s)] is 
obtainable: 

rc(s)] = [G(s)][P(s)][E(s)] (2.1) 

where the error or actuating function [E(s) ] yields 

[E(s)I - [R(s)j - [H(s)][C(s)] ; (2.2) 

the notation is defined in the S5rmbols table. Substituting equation 
(2.2) into (2.1) and solving for [C(s)3 yields the output relation, 

rc(s)j = [G(s)][P(s)]{H- [H(s)][G(s)][P(s)]r'[R(s)I (2.3) 

Equation (2.3) gives the closed-loop- output response in terms of the 
input vector. Suppose that the kth diagonal element of [H(s) ] is 
set to zero and all the elements of [R(s)] are nulled except the kth 
element which is set to unity; the result is the frequency response 
between the kth input and the outputs when the kth loop is open. To 
simplify notation, define 

[V(s)3 4 {I + [H(s)][G(s)][P(s)]}"' (2.4) 

and 

[U(s)l A [G.(s)][P(s)I . (2.5) 

Hence the open-loop complex frequency response of the kth system is 

== u^(s) • 


Cj^(s) 

H,(s) 


( 2 . 6 ) 
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R(s) IS THE N X 1 INPUT VECTOR 

E(-s) IS THE N X 1 ERROR VECTOR 

C(s) IS THE N X 1 OUTPUT VECTOR 

P(s) IS THE M X N MATRIX TRANSFER 

FUNCTION DESCRIBING THE PLANT 
G(s) 'IS THE N X M MATRIX TRANSFER 

FUNCTION DESCRIBING THE CONTROL 
LAW AND CASCADED COMPENSATION 
H(s). IS AN N X N UNITY MATRIX 


Figure 4. Vector Representation of the 
Multivariable Control System. 
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where la (s) and Vj^(s) are, respectively, the kth row of [U(s)] and 
the kth column of [V(s)] with the kth diagonal element of [H(s)] 
set to zero. By fixing the proper diagonal element of tH(s)] to 
zero, it is obvious how equations (2,4), (2.5), and (2.6) can be 
used to elvaluate the open-loop frequency response for each kth 
system. 


Evaluation of the Critical frequencies 

Next attention is focused on the determination of the critical 
frequencies with respect to the design objectives. In CIP the 
design is accomplished by requiring the open-loop frequency .response 
to satisfy certain specifications. These design specifications are 
converted to distances between certain critical points of the open- 
loop frequency response and certain points in the corresponding 
complex plane where s = jw. The typical objective function for the 
kth open-loop system is 

d = j [A + Cj^(ju)3[A + Cj^(jw)]*|^ (2.7) 

where A is the point in the complex plane from which the specifica- 
tion is measured; e.g. , for stability margins A is the (-1 + jO) 
point. In equation (2.7), the response Cj^(jw) is calculated from 
(2.6) with Rj^(ju) set equal to unity. 

Calculation of the Partial Vectors 

Now attention is directed to the calculation of the change in 
the design objectives with respect to the free parameters of the 
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controller. With this objective in mind, the partial derivative of 
the distance d with respect to some parameter w is 


Sd 1 \ 3cJ(jO)) ac (jo)) J 

^ - 2 [A+ 


^d . 


Equation (2.8) can be rewritten as 
If - Re 


3Cj^(jw)) 


aw ^ 


^d 


( 2 . 8 ) 


(2.9) 


Evaluation of (2.9) depends on determining accurately the partial 
term, 3Cj^(jO)) /aw. Using the chain rule this becomes 


aCj^(ju) 9Cj^(ju) ac^^ 

3w 8G. . 9w 

11 


( 2 . 10 ) 


where G. . is the element of the controller [G(s)] in which the 
ij 

free parameter w appears with i = 1,2,...,N system outputs and 
j = 1,2,...,M system states sensed. 

The necessary equations for evaluating the first term in (2.10) 
are derived in the sequel. The partial of (2.3) with respect to 
the controller element G. . gives 


= ~[G(s)][P(s)]{l + [H(s)][G(s)][P(s)]}"^[H(s)l • 

[P(s)3{l + [H(s)][G(s)][P(s)]p[R(s)3 + 


[P(s)l{l + [H(s)][G(s)][P(s)l}”^[R(s)] . (2.11) 


Then, the 3Cj^(s)/aG^, is the kth element of (2.11) with the kth 
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diagonal element of [H(s) ] set to zero, and all the elements of 
[R(s)] set to zero except the kth which is set to unity. In (2.11) 
the partial term B[G(s) ]/9G^^ is 8(zero matrix except for the (ij)th 
element which is unity. 

Next consideration is given to the evaluation of the second 
term in equation (2.10). It is assumed that the (ij)th element of 
[G(s) ) is composed of a cascaded arrangement of transfer functions, 

i.e. , 

K 

G (s) = n G..,(s) (2.12) 

where K is the number of cascaded elements. The £th cascaded ele- 
ment of the (ij)th element of [G(s)] has the general form 


M 


■ ¥ 




m 


(2.13) 


Jo 


n 


Then, the free parameters of this element are the x’s and y’s. If 
w in (2.10) is the pth numerator coefficient of (2.13), then 


or 


3G^.(s) 
9x , » 

ij£p 


+ sP 

N 


r n 

io 


9G. .(s) 

--iJ 

9x. 

XjSrp 


G„(s) 


+sl’ 


M 

I 

m=0 


n 


X. , „ s 
xjS,m 


(2.14) 


(2.15) 
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Similarly if w in (2.10) is the pth denominator coefficient of 
(2.13), then 

3G..(s) K 

— = n 

^^ij£p k=l 

or 

9G..(s) 

U 

By appropriately using equations (2.3), (2.9), (2.10), (2.12), 
(2.15) and (2.17), the first order change of any CIP objective 
function with respect to the free parameters of the controller can 


~S‘ 


'^ij^^^ / N 


(2.17) 


V 

n=0 


ijk 


(s) 


-s" 


N 


i.n-0 


(2.16) 


^ijiln" 


n 


be calculated. 



CHAPTER III 


DESCRIPTION OF THE COMPUTATIONAL FLOW DIAGRAM 
Implementation of the Algorithm 

The algorithm implementation evolved with the following 
objectives in mind. 

1. Any linear, time-invariant system structure should be 
acceptable. 

2. Numerical problems' should not restrict the method to low 
order systems. 

3. Computational requirements should be reasonable for high 
order systems . 

4. Monitoring of compensation at desired iteration levels 
should be possible. 

■ 5. Partial derivatives should be exact. 

A Synopsis of the Algorithm 

Keying on the aforementioned goals , possibilities were weighed 
to determine the most effective methods of implementing and computer 
coding the algorithm. The following is a description of the princi- 
ples and computational logic involved in producing an executable 
version of the compensator improvement program for the multivariable 
control system of Figure 2. In this sequel the- logic flow diagram 
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of the Continuance Criterion, (D) The Determination of the Gradient 
Vectors Corresponding to the Active Constraints, and (E) Calculation 
of the Directional Vector and the Compensation Enhancement. For a 
detailed explanation of the particular subroutines used, refer to 
alphabetical listing in Appendix A. The Fortran IV computer code is 
listed alphabetically in Appendix B. 


Section A; Data Input 

Referring to the schematic diagram of Figure 2, recall that the 
multivariable system may be viewed as n coupled feedback systems — 
one for each element of the input vector. With this view, the loop 
transfer function of the kth system is obtainable by opening the 
feedback path at and determining the response Cj^(s)/R^(s) with 
all input R's set to zero except the kth element. Thus the variable 
k in the computational flow diagram of Figure 3 is defined as the 
respective subsystem in accordance with the corresponding element of 
the input vector. 

The input description of the multivariable configuration 
requires discrete frequency data from each input to each output in 
describing the plant system; whereas, the initial compensation for 
each controller is described by a matrix of transfer functions. 

With this information the open-loop frequency response is obtained 
for each of the n coupled systems. Likewise, for each subsystem, a 
set of critical points, that is, frequencies at which margins of 
stability or attenuation occur, is determined. Hence, the input 
routine requires data of four types as shown in Table 1 and clari- 
fied in the following discussion. 



1. Iteration Control 

a. Mode, identification code 

b. Start, stop, print iterations 

c. Maxiroum, ■minimum step sizes 
Etc. 

2. Design Specifications 

a. Desired Stability and Attenuation Radii 

b. Frequency Ranges over which searches for 
critical points are to be made 

3. Description of Plant 

a. Number of Control Inputs 

b. Number of Outputs 

c. Discrete frequency response data 

4. Description of Compensation 

a. Gain Constant in each channel 

b. Number of Subcompensators in each channel 

c. Coefficients for each subcompensator in first 
and second order factors only 

d. Constraints to be placed on the coefficients 
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1. Iteration Control 

First, user control parameters are entered; these include the 
extremum step sizes to be taken on iterations, maximum iterations 
for convergence, designation of iterations to be printed, user 
identification code, etc. Here also the user must specify the mode 
used in the program to determine when an iteration has been com- 
pleted, In particular, the mode designates which continuance 
criterion must be used to determine whether the trial design at the 
Ci + l)th iteration is an improvement in comparison to the results 
at the ith iteration. One of two modes must be chosen; 

i. Total Improved Frequency Response Mode (TIER) re- 
quires that from iteration to iteration no unsatis- 
fied objectives or design specifications are allowed 
to degrade and insures improvement in at least one. 
il. Sum Improved Frequency Response Mode (SIFR) requires 
that the sum improvement exceed the sum degradation 
from iteration to iteration. 

It is obvious that the TIFR mode produces a more stringent continuance 
criterion on the compensation. 

2. Design Specifications 

The second portion of the input data 'designates the design 
specifications for achieving relative stability and relative 
attenuation. In particular, the mathematical formulation of the 
design problem is to determine the free parameters of the compen- 
sators such that the objective functions satisfy a set of design 
specifications, i.e., gain, phase, stability and attenuation 
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margins. Mathematically, if • a total of ti critical frequency response 
points are chosen, the problem can be expressed as a strict con- 
straint mathematical programming problem of the form; 

T 

Determine the 3c 
T 

such that the constraints g.(x ) > b. , i == l,2,,,,,n : (3-sl) 

X “ — X 

where x represents the free compensator parameters, b^ represents 

T 

the design specifications, and g (x ) contains the objective func- 

X ~ 

tions , that is , frequency response limitations and constraints . 

Thus the general idea is to change the compensator coefficients so 
that each constraint comes closer to being satisfied at each itera- 
tion.' Other methods of obtaining the design objectives could be 
implemented; however, from a practical point of view the method of 
the strict constraint problem is particularly appealing in pro- 
ducing a change vector for the compensator coefficients that 
insures the existence of a Total and/or Sum Improved Response. 
Furthermore, this method allows the margin radii specifications to 
become frequency dependent. Conceivably, it is desirable that 
regions of the frequency response be various distances from the 
(-1 + jO) point in the GH(jto) plane while other regions be con- 
strained to be greater or less than limitations with respect to the 
origin of the GH(joi) plane. Thus in general a frequency response 
is desired to have some basic shape which can be translated with 
respect to frequency and not constrained to match exactly a desired, 
frequency response. 
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3. Description of Plant 

The multivariable configuration of Figure 2 may be described 
in two parts, the plant and the controller. First the description 
of the uncompensated plant requires discrete frequency response 
data between each input and each output channel. The choice of the 
discrete data description for the plant was made to avoid compu- 
tational difficulties that might be encountered in evaluating high 
order transfer functions and to conserve computing time in the 
iterative process of CIP. Furthermore, discrete frequency response 
data is often the best information available for describing the 
system. 


4. Description of Compensation 

Secondly, the initial compensation or controller design is 
necessary for each control input. Again recall the objective of 
this work is not to develop a self-contained, computer-aided design 
algorithm but to provide the control engineer with a design aid. 

In this regard, it is assumed, that the designer knows the control 
law necessary to enhance his design objectives. Normally the engi- 
neer uses s-domain rational functions in investigating designs, 
thus for simplicity, the compensation elements are described by 
transfer, functions in the form of cascaded first and second order 
factors, that is,. 


G(s) 


(GAIN) 


N1 • N-2 

■n (ZA, + ZB.s) n (ZC. + ZD.s + ZE.s^) 
i=l ^ ^ j =1 ^ j J 

Ml M2 ' 

H (PA. + PB.s) n (PC. + PD.s + PE.s^) 
i=l ^ ^ j=l J J J 


(3.2) 
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Section B; Frequency Response Manipulations 

Figure 5 gives a detailed expansion of the logic of Section B 
in Figure 3. Note that Section B is composed of four major rou- 
tines in determining the frequency response: Delete Points, Add 

Points, Calculation of the Closed-Loop Frequency Response, Calcu- 
lation of Critical Points. 

The Delete Points routine is designed to remove any frequency 
points and corresponding response terms no longer of major concern 
which might have been added for accuracy on previous iterations; 
however, the original data are always retained. The routine is 
designed to save computer storage as well as' computer time in 
response calculations and in scanning for margins. ■ This algorithm 
is coded in the subprogram DELETE [10]. 

The Calculation of the Closed-Loop Frequency Response is per- 
formed by implementing the equations developed in Section 3 of 
Chapter II. In particular, the closed-loop output vector [C(s)] is 
determined by the relation, 

[C(s)l = [G(s)][P(s)]{l + [H(s)][G(s)][P(s)]} -l[R(s)] . (3.3) 

Equation (3.3) gives the closed-loop output vector in terms of the 
input vector. ' Based on the theoretical concepts of Chapter II, 
suppose that the kth diagonal element of [H(s)] is set to zero and 
all the elements of [R(s) ] are nulled except the kth element which is 
set to unity; the result is the frequency response between the kth 
input and the outputs when the kth loop is open. Utilizing the 


aforementioned notation, 



A 

DATA 

INPUT 


SET K = 0 



Figure 5. Logic Diagram Representing the Frequency 
Response Manipulations of Section B. 
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[V(s)] A {l + [H(s)][G(s)][P(s)]}"' (3.4) 

and 

tU(s)l A [G(s)][P(s)] . (3.5) 

The open-loop frequency response of the kth system is 

= £^(s)vj^(s) , (3.6) 

where ja (s) and are, respectively, the kth row of [U(s)], 

and kth column of [V(s)]. By setting the proper diagonal element 
of [H(s) ] to zero, equations (3.4), (3.5), and (3.6) can be used to 
evaluate the frequency response of each system. 

This algorithm is coded in the main program using two sub- 
routines: EVAL and CRT. The subprogram ,EVAL evaluates the con- 

troller at the specified frequencies with the aid of the program 
POIEV, a polynomial evaluation routine, EVAL then detemines the 
product of the controller response and the plant response, that is, 
[G(s)l • [P(s)], With this transfer relation the subprogram CRT 
determines the total response [C(s)] and selects the open-loop 
frequency response of the kth system. 

The Determination of the Critical Points is designed to yield 
the critical points of the open— loop frequency response and to 
ascertain whether they satisfy certain design specifications for 
achieving the relative stability and relative attenuation margins . 
Recall these design specifications are expressed mathematically as 
the distances between certain critical points of the frequency 


R^(s) 
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response and particular points in the corresponding complex plane. 
The typical objective function for the kth open-loop system is 

d = |[A + C^(joj)][A + C^Cjco)]*!'^ (3.7) 

where A is the point in .the complex plane from which the specifica- 
tion is measured; for stability, gain, and phase margins A is the 
(-1 + jO) point; for attenuation margins, the (0 + jO) point is 
chosen. The design specifications include subprograms for deter- 
mining the gain, phase, stability, and attenuation margins. 

The subprogram .to Add Points is designed to provide more data 
around each of the critical frequency points, thereby, yielding a 
more exact margin value without the input of excessive data and the 
consequent increase in storage. This algorithm is encoded as ADOPTS 
and uses an interpolate routine INTER, as well as the aforementioned 
EVAL and CRT routines to update the frequency response for each 
added data point. ■ 

As indicated in Figure 5, the frequency response and critical 
margins are calculated for each kth system adding and deleting 
frequency points as necessary. 

Section C: Evaluation of the Continuance Criterion 

Figure 6 represents the logic decision blocks of Section C in 
Figure 3. These decision blocks are encoded in the main program 
and are designed to force the program into the specified continuance 
criterion when an iteration has been completed. Recall that the 
continuance criterion mode must be specified by the user as input 




3 
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data to determine whether the trial design at the (i + l)th 
iteration is an improvement in comparison to the results of the itli 
iteration. The two modes available, the TIFR mode and the SIFR mode, 
are as defined in Section A. 

Note that for the first iteration the mode block is bypassed 
thereby assuming that an improved solution has occurred and allowing 
the program to continue to the determination of the partial vectors. 

If the SIFR/TIFR condition for the current iteration is satis- 
fied the program checks user control data to decide whether the 
maximum iteration condition has been exceeded. If the last 
iteration has been reached the program sets a stop condition which 
prohibits further manipulations of the active constraints and partial 
vectors. Assuming, the maximum iteration, code has not been met, the 
main program directs control to Section D. 

Now if the SIFR/TIFR criterion has not been met, the program 
interprets this condition to mean that the change in the compensator 
coefficients was too large and control proceeds to decrease the step 
size of the previous iteration in an attempt to force an improved 
solution. 

Section D; Determination of the Gradient Vectors 

Figure 7 -is an expanded view of Block D of Figure 3; Section D 
,is concerned with the determination of the active constraints, that 
is, the. margins which do not satisfy the required design specifica- 
tions, and their relation in evaluating the partial vectors. 

The selection of active constraints is coded within the main 
program, CIP checks to determine which of the specified stability 
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and attenuation margins do not satisfy the margin design specifi- 
cations entered by the user as input data. As previously noted, 
it is possible to demand any combination of margin requirements as 
gain, phase, stability, and attenuation radii. These margins can 
be manipulated so as to make the specifications or constraints 
frequency dependent. A list of the, active radii requirements, 
margins, and corresponding frequencies is then prepared to alleviate 
any margins already satisfied. The program checks for any dupli- 
cations in margins, in which case retaining only the first such 
critical point found. 

The second objective of Section D is the calculation of the 
partial vectors as described in Section 3 of Chapter II. Recall 
that the partial vectors represent the change in the design objec- 
tives with respect to the free parameters, that is, the compensator 
coefficients of the controller. Thus the partial derivative of the 
objective function d with respect to some parameter w can be 
expressed as 

c, (jw) \ 

II =.Re|[A-f C^(jU))]*-^^} ^ d . (3.8) 

Accordingly, the partial term 9C (ju)/3w can be expanded by 
the chain rule as 


aCj^Cjw) aCj^Cjw) 3G_ 
3w 3G . . 3w 


(3.9) 


where G.,. is the element of the controller [G(s)] in which the free 


parameter w appears. 
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Evaluation of the first term in (3. 9) yields the relation 


9[C(s)] 

9Gy 


- [G(s)][P(s)]{l + [H(s)][G(s)][P(s)]}"^tH(s)] • 

[P(s)]{l + tH(s)ltG(s)][P(s)]r'tR(s)l + 
ij 


[P(s)l{l + [H(s)][G(s)][P(s)]r'[R(s)] (3.10) 

ij 


Then, 3Cj^(s)/9G^^ is the kth element of (3.10) with the kth 
diagonal element of lH(s)] set to zero, and all the elements of 
[R(s)l set to zero except the kth which is set to unity. In (3.10) 
3 [G(s)]/ 9G^^ is a zero matrix except for the (ij)th element which is 
unity. 

Evaluation of the second partial term in equation (3.9) is 
acquired by assuming that the (ij)th element of the compensator 
matrix [G(s)] is composed of a cascaded-arrangement of transfer 
functions, i.e., 

K 

G (s) = n G (s) (3.11) 

11 . IJR 


where K represents the number of cascaded elements. Thus the £th 
cascaded element of the (ij)th compensator has the general form 


M 


I 


m 




m=0 


ijAm" 


io 


(3.12) 


where the free parameters of this element are the x’s and y’s. 
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Assuming the parameter w in (3.9) is the pth numerator coef- 
ficient in (3.12), then the following expression is obtained: 


3G..(s) 
. ^3 , , 
3x. 

1J&P 


+ s' 




(3.13) 


/ X. . „ s 
I*- xjAra 
m=0 


Similarly, letting w represent the pth denominator coefficient 
in equation (3.9), then 

9G. . (s) 

3-3 " 

3y - * 0 

■"xjS-p 

Thus equations (3.10), (3.13), and (3.14) can be used effective- 
ly to determine the, first order change of any GIF objective function 
with respect to the free parameters of the controller. 

As indicated by the decision block, the partial vector of each 
system is tabulated and stored in an orderly array for later use in 
determining the directional vector in Section E. The partial vector 
routine is coded in the subprogram PARTAL in conjunction with the 
frequency response subroutine CRT; subroutine CRT determines the 

partial term 3C, (ja))/3G, . in equation (3.9). 

K- xj 

Section E; Calculation of the Directional Vector and 
Compensation Enhancement 

Of major significance in Section E is the manipulation of the 
system partial vectors in obtaining a directional vector using the 
constraint improvement algorithm. (See Figure 8.) 


^ij^^^ / N 


(3.14) 


V n ’ 

h ^i3an" 

n=0 I 
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Figure 8. Logic Diagram Representing the Calculation of the 
Directional Vector and Corresponding Compensator 
Enhancement of Section E. 
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First the decision block of Section E tests to ascertain 
whether the stop condition has been set. If so, the output routine 
OUTPT is called; otherwise, the main program determines whether the 
design specifications have been met and calls the output routine if 
necessary. Assuming the design objectives are not satisfied, the 
program checks user data to determine whether an output is desired 
for the particular iteration. The directional vector then is cal- 
culated with the aid of the Constraint Improvement Algorithm [7] in 
the subprogram pIRVEC. This subprogram also checks for routine 
failure yielding a stop command. The subprograms MATMUL and MATINV 
are auxiliary routines to DIRVEG for matrix multiplication and 
inversion, respectively. 

In Figure 8, after calculating the directional vector and 
assuming the CIA did not fail, CIP investigates the user* option of 
constraining the poles and zeros of compensation to lie within a 
specified damping ratio sector. By limiting the compensation to 
first and second order factors, the complexity of constraining the 
compensation poles and zeros to' lie within a sector defined by con- 
stant damping ratio lines as shown in Figure 9 is reduced greatly. 
Hence the next step is the determination of the directions of move- 
ment of the compensator poles and zeros on the specified zeta 
boundaries. In order to avoid a zeta violation, the directions of 
these poles and zeros must.be along the boundaries or into the 
defined sector. If with the present directional vector the move- 
ments of these poles and zeros are in the wrong directions, judi- 
ciously selected terms in the partial vectors are nulled and the 



Allowable 
Zeta Sector 



Figure 9 . s~Plane Configuratior 2Ctor 

of Desirable Gompenss 3 and 

Typical Pole Locationc j.wi. a der 

System. 
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directional vector is recomputed. This process is continued until 
the directions of movement of all compensation poles and zeros on 
the boundaries do not result in a zeta violation. The checking of 
the directions of movement of these poles and zeros and the setting 
of the elements of the partial vectors to zero is accomplished by 
the subprogram XCHECK. This routine assures the existence of a non- 
zero step size that will not produce a zeta violation by poles and 
zeros on the boundaries. 

Referring to Figure 8, after an acceptable directional vector 
has been established in accordance with XCHECK, a step size is 
selected, and in conjunction with this directional vector, the 
individual compensator coefficients are effectively augmented. At 
this point, a violation in the zdta constraints can occur from an 
inappropriate selection of the step size, i.e. , usually if the step 
size is too large. Thus, the zeta constraints are checked. If a 
violation occurs, the maximum step size that will not produce a 
violation is computed and the compensator coefficients are rein- 
cremented; otherwise, the program recycles to Section B in Figure 3. 
The check for zeta violations, as well as the computation of a 
maximum acceptable step size, is accomplished by the subroutine 
YCHECK. The theory underlying this routine along with additional 
usage information is presented in Appendix A, 



CHAPTER IV 


INVESTIGATIONS AND EXAMPLES 

In order to illustrate the practical utility of the multivaria- 
ble Compensator Improvement Program, the Improvements of the com- 
pensators for space related examples are presented. This by no means 
limits the scope of the work to space oriented control systems, but 
rather provides large system problems which have been investigated 
by other means. In particular, three examples are discussed: (1) a 

dual input, dual output system with uncoupled characteristics; 

(2) a dual input, dual output system with coupling; (3) a dual input, 
four output system exhibiting coupled characteristics. 

Uncoupled Dual Input, Dual Output System 

In this example the system under consideration is similar to 
that of Figure 2 with M=2 controller inputs or measured states, and 
N=2 controller outputs. Figure 10 shows the actual subsystem under 
investigation and is representative of the attitude control system 
for a finned launch vehicle at a specified flight time following 
launch [11]. Each subsystem has plant dynamics C^(s)/R^(s). 
described by the uncompensated frequency response plot of Figure 11 
where k represents the number of controller inputs and hence the 
number of subsystems; k equals two in this example. Table 2 
exhibits the twenty-eight discrete frequency points chosen to 
describe the open-loop response of each subsystem. The compensation 
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Figure 10. Block Diagram Representing Each Subsystem 
of the Finned Vehicle Example. 


Table 2. Frequency Response Input Data Describing 
the Plant in the Finned Vehicle Example. 


Data 

Points 

Complex Frequency 
:RE[s] IM[s] 

Uncompensated Plant Response 
RE[P] • IM[P] 

1 

0.000 

0.100 

-69,6000 

54.9000 

2 

0.000 

0.132 

-67,3000 

36.6000 

3 

0.000 

0.178 

-63.3000 

20.8000 

4 

0.000 

0.240 

-57,2000 

8.9000 

5 

0.000 

0.347 

-46.9600 

- 0.5890 

6 

0.000 

0.501 

-34.8000 

- 4.3200 

7 

0.000 

0.646 

-26.6000 

- 4.3300 

8 

0.000 

1.230 

-10.8800 

- 1.5220 

9 

0.000 

, 2.000 

- 4.7100 

- 0.3060 

10 

0.100 

2.800 

- 2.4800 

- 0.1840> 

11 

0.200 

3.600 

- 1.5100 

- 0.1020 

12 

0.250 

4.000 - 

- 1.2200 

- U.0750 

13 

0.250 

4.921 

- 0.8079 

0.0018 

14 

0.250 

6.000 

- 0.5369 

0.0340 

15 

0.200 

6.400 

- 0.4690 

0.0469 

16 

•0.100 

7.200 

- 0.3630 

0.0603 

17 

0.000 

8.000 

- 0.2870 

0.0645 

18 

0,000 

10.071 

- 0.1666 

0.0478 

19 

0,000 

12.400 

- 0.0956 

0.0309 

20 

0.000 

15,600 

- 0.0402 

0.0013 

21 

0.000 

18.327 

- 0.5180 

- 0.1090 

22 

0.000 

19.191 

- 0.1924 

- 0.0525 

23 

0.000 

19.638 

- 0.1677 

0.0249 

24 

0.000 

20.563 

- 0.0969 

0.0534 

25 

0.000 

30.415 

- 0.0139 

0.0182 

26 

0.000 

40.095 

- 0.0047 

0.0092 

27 

0.000 

60.686 

- 0.0006 

0.0031 

28 

0.000 

100.000 

- 0.0001 

0.0007 
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matrix consists of identical compensators in tlie diagonal elements, 

•L. • j 


Gji(s) - 622(3) 


(1 + 10s) (5 + 1.25s) 
(1 + 11s) (5 + 1.00s) 


(4.1) • 


whereas, the off-diagonal terms are chosen as zero to inhibit any 
coupling terms. 

It is desired to modify the compensators, G^j(s) and 6^2 (s)» 
so that the closed-loop step response of each subsystem reasonably 
is damped and "ringing" caused by the low-damped high frequency modes 
is negligible. Further, the DC gain of each compensator is chosen 
so that the magnitudes of the steady-state errors to a velocity 
input are less than 0.15. These specifications are satisfactorily 
achieved by requiring that 


(1) 

all SM' s ^ 

0 . 5 when 0 


hi 

< 16.0 

(2) 

all AM's _< 

0.1 when 16 

< 

W 

£ 100 

(3) 

the DC- gain 

of the compensator 

is 

greater 


than 26.67. 






After 31 iterations, approximately 36 seconds of CPU 
Univac 1108 Computer, the design compensation is obtained 


time on a 
as 


Gii(s) 


622(3) 


(1.0 + 1.60084s) (5.0 + 5.57314s) 
(1.0 + 16.3982s) (5.0 + 1.14051s) 


(4.3) 


The compensated frequency response for each kth subsystem is 
illustrated in Figure 12 in which all design specifications have 
been accomplished. In particular. Entry A in Table 3 shows the 
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Figure 11. Frequency Response Representing Each Susbsystem of the 
Uncompensated System of the Finned Vehicle Example. 


90 * 



270 * 

Figure 12. The Improved Compensated Frequency Response of Example 1. 
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Table 3. System Specifications for the Finned Vehicle Examples. 


Margin 

Number 

Margin 

Value 

Complex Frequency 
RE[s] IMts] 

Desired 

Margin 

Margin 

T)rpe 

Active 

List 



A. The Uncoupled 

System: 



Iteration No. 0 


Subsystem 1, 

2. 



1 

42.2900 

0.000 

0.347 

0.50 

GM 

No 

2 

0.5036 

0.200 

6.400 

0.50 

GM 

Yes 

3 

7.5340 

0.250 

4.437 

30.00 

PM 

Yes ' 

4 

0.2241 

0.000 

19.190 

0.10 

AM 

Yes 

Iteration No. 31 






1 

3.1130 

. 0.000 

0.673 

0.50 

GM 

No 

2 

30.0300 

0.000 

2.086 

30.00 

PM 

No 

3 

■ 0.0920 

0.000 

19.190 

0.10 

AM 

No 



B. The Coupled System: 



Iteration No. 0 


Subsystem 1. 




1 

42.2900 

0.000 

0.347 

0.50 

GM 

No 

2 

0.6114 

0.000 

8.000 

0.50 

GM 

Yes 

3 

11.0500 

0.250 

4.613 

30.00 

PM 

Yes 

4 

0.2976 

0.000 

19.190 

0.10 

AM 

Yes 




Subsystem 2, 




1 

29.3000 

0.000 

0.347 

0.50 

GM 

No 

2 

0.2010 

0.200 

3.600 

30.00 

PM 

Yes 

3 

0.2360 

0.000' 

19.190 

0.10 

AM 

Yes 

Iteration No. 50 

Subsystem 1. 




1 

2.4590 

6.000 

0.700 

0.50 

GM 

No 

2 

28.0300 

0.000 

1.882 

30.00 

PM 

No 

3 

0.0870 

0.000 

19.190 

0.10 

AM 

No 



Subsystem 2. 




1 

30.0880 

0.000 

0.606 

0.50 

GM 

No 

2 

29.4700 

0.000 

1.882 

30.00 

PM 

No 

3 

0.0960 

0.000 

19.190 

0.10 

AM 

No 
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specified design objectives versus the final system specifications. 
Similar results were obtained with the Compensator Improvement 
Program by Mitchell and McDaniel in reference [11] for the single 
input, output system. 


A Coupled Dual Input, Dual Output System 


This example utilizes the same system as the previous illus- 
tration but in this case the plant subsystems are coupled with non- 
zero off-diagonal terms. Thus, the same uncompensated frequency 
response data is used to describe the diagonal terms of the plant 
matrix [P(s)]*. Table 4 gives the frequency, response data des- 
cribing the coupling terms, that is, compen- 

sation matrix remains the same as described in equation (4.3); 
however, for generality, the compensator gain of the G^^(s) element 
has been altered to a factor of 0.7 instead of unity. 

Requiring the same design objectives as stated in equations 
(4.2), the improved compensation elements, 


(1.0 + 1.42713s) (5.0 + 6.52626s) 
(1.0 + 19.2239s) (5.0 + 1.25501s) 


(4.4) 


and 


^22 (s) 


0. 7(1.0 + 2.35364s) (5.0 + 4.07685s) 
(1.0 + 15.2297s) (5.0 + 1.01990s) 


(4.5) 


are obtained after 50 iterations. Entry B in Table 3 shows the 
desired objectives as compared to the final system specifications 
after the compensation improvement. 
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Table 4. Frequency Response Data Describing the Coupling 

Elements of the Plant Matrix [P(s)] in the Coupled 
Finned Vehicle Example 


Data Complex Frequency Uncompensated Plant Response 


Point RE[s] IM[s] RE[P^^3 IM[P^^] RE[P^^] 


1 

0.000 

0.100 

0.1000 

0.0005 

0.2001 

0.0066 

2 

0.000 

0.132 

0.1000 

0.0066 

0.2002 

0.0088 

3 

0.000 

0.178 

0.1000 

0.0009 

0.2003 

0.0118 

4 

0.000 

0.240 

0.1001 

0.0012 

0.2006 

0.0159 

5 

0.000 

0.347 

0.1001 

0.0017 

0.2013 

0.0230 

6 

0.000 

0.501 

0.1003 

0.0029 

0.2028 

0.0331 

7 

0.000 

0.646 

0.1004 ■ 

0.0032 

0.2046 

0.0425 

8 

0.000 

1.230 

0.1014 

0.0058 

0.2161 

0.0786 

9 

0.000 

2.000 

0.1035 

0.0086 

0.2400 

0.1200 

10 

0.100 

2.800 

0.1062 

0.0134 

0.2750 

0.1492 1 

11 

0.200 

3.600 

0.1087 

0.0112 

0.3105 

0.1681 

12 

0.250 

4.000 

0.1099 

0.0115 

0.3276 

0.1744 

13 

0.250 

4.921 

0.1123 

0.0119 

0.3630 

0.1866 

14 

0.250 

6.000 

0.1147 

0.0118 

0.4002 

0.1988 

15 

0.200 

6.400 

0.1154 

0.0118 

0.1259 

0.1935 

16 

0.100 

7.200 

0.1168 

0.0115 

0.4356 

0.1940 

17 

0.000 

8.000 

0.1176 

0.0102 

0.4520 

0.1742 

18 

0.000 

10.071 

0.1201 

0.0099 

0.4952 

0.1759 

19 

0.000 

12.400 

0.1215 

0.0087 

0.5241 

0.1568 

20 

0.000 

15.600 

0.1227 

0.0073 

0.5485 

0.1340 

21 

0.000 

18.327 

0.1233 

0.0063 

0.5613 

0.1183 

22 

0.000 

19.191 

0.1234 

0.0061 

0.5644 

0.1139 

23 

0.000 

19.638 

0.1235 

0.0059 

0.5658 

0.1118 

24 

0.000 

20.563 

0.1236 

0.0057 

0.5686 

0.1076 ; 

25 

0.000 

30.415 

0.1243 

0.0040 

0.5850 

0.0759 

26 

0.000 

40.095 

0.1246 

0.0030 

0.5912 

0.0585 

27 

0.000 

60.686 

0.1248 

0.0020 

0.5061 

0.0391 

28 

0.000 

100.000 

0.1249 

0.0012 

0.5986 

0.0239 
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A Dual Inputs Four Output System with 
Coupled Characteristics 

This example is representative of the Yaw/Roll Ascent Flight 
Control System for the Space Shuttle. The system is similar to that 
of Figure 2 with a plant possessing two control inputs and four out- 
puts. The 23 discrete frequency response data chosen to describe 
the plant dynamics are listed in Table 5. The compensation matrix 
[G(s)] given in Table 6 actually is designed for use on the space 
shuttle and exemplifies the complexity required in achieving a set 
of design objectives. 

Given the design requirements of Table 7, the CIP produced the 
improved compensation matrix of Table 6 in 5 iterations, that is, 

20 seconds of CPU time on the Univac 1108. 

In summary, the CIP is a fast and effective design tool in the 
area of compensation improvement. For the Finned Vehicle examples, 
approximately 28K words of core storage is required; the Shuttle 
example executes in 32K of storage. In each example the programming 
time and computer time is minimal. 
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Table 5. Plant Dynamics Describing the Yaw/Roli 
Ascent Flight Control System for Space 
Shuttle Example. 


Data 

Frequency 




Uncompensated Plant Frequency Response Data 



Point 

RE[s] 

RE[Pj^l 


RElP^il 

RE[P22l 

RE[P33^) 

REtP32l 

RE[P^il 



IM[s] 


IM[P^2l 

IM[p2i] 

IMtp22‘] 

im[P3i1 

IM[P32l 

IM[P,ll 

IM[P„] 

1 

0.0000 


2.8340 

- 0.0431 

- 10.6900 

1.0700 

2.7380 

0.1860 

31,2900 

- 5.8820 


0.0100 

- 

3.2070 

- 0.0463 

- 2,3320 

- 0.1385 

-13.7400 

2.1720 

22.7100 

- 2.7820 

2 

0.0000 

“ 

0.7722 

- 0.0117 

- 2.6520 

1.1850 

- 0.1172 

0.1527 

3.6560 

- 1.9500 


0.0185 

- 

1.1690 

- 0.0171 

- 0.1011 

- 0.2079 

- 2.6890 

1.2350 

1.9790 

- 1,6310 

3 

0.0000 

- 

0.6223 

- 0.0095 

- 2.1060 

- 1.1910 

- 0.1571 

0.1525 

2.5500 

- 1.5740 


0.0209 

- 

1.0030 

- 0.0147 

- 0.1604 

- 0.2360 

- 2.0880 

1.0930 

1.3180 

- 1.4430 

It 

0.0000 

- 

0.4073 

- 0.0063 

- 1.3360 

1.1950 

- 0.1815 

0.1522 

1.2320 

- 1.0060 


0.0275 

- 

0.7280 

- 0.0108 

0.2186 

- 0.3130 

- 1.2750 

0,8303 

0.5690 

- 1.0950 

5 

0.0000 

- 

0.1988 

- 0.0036 

- 0,6279 

1.1670 

- 0.1668 

0.1486 

0.3088 

- 0.4414 


0.0507 

- 

-.3513 

- 0.0059 

0.2656 

- 0.5978 

- 0.5097 

0.4312 

0.1152 

- 0.5680 

6 

0.0000 

- 

0.1948 

- 0.0045 

- 0.S711 

1.1680 

- 0.1802 

0.1479 

0.3197 

- 0.4330 


0.0517 

- 

0.3422 

- 0.0059 

- 0,0029 

. 0.6239 

- 0.4985 

0.4219 

0,1078 

- 0.5562 

7 

0.0000 

- 

0.1940 

- 0.0049 

- 0.8675 

1.1520 

- 0.1878 

0.1476 

0.3208 

- 0.4318 


0.0519 

- 

0.3403 

- 0.0049 

- 0.1351 

- 0.6331 

- 0.4795 

0.-4205 

0.0795 

- 0.5532 

8 

0.0000 

- 

0,6462 

- 0.2679 

- 1.2820 

1.1400 

0.6761 

0.1210 

- 0.0047 

- 0.3364 


0.0659 

- 

0.1136 

- 0.0085 

5.0450 

- 0,9269 

- 0.0904 

0.3103 

- 0.3594 

- 0.4264 

9 

0.0000 


0.5705 

- 0.0246 

4.4650 

0.9682 

0.7388 

0,1190 

0.4012 

- 0.3467 


0.0662 

- 

1.1130 

0.0212 

4.4600 

- 0.9142 

- 1.1540 

0.3393 

0.4899 

- 0.4268 

10 

0.0000 

- 

0.1037 

- 0.0046 

3.7080 

0.9895 

0.0113 

0.1401 

0.4222 

- 0.3450 


0.0665 

- 

0.9970 

0.0179 

0,5986 

- 0.8032 

- 1.1280 

0.3372 

0.1397 

- 0.4138 

11 

0.0000 

- 

0.2463 

- 0.0002 

0.4308 

1.0820 

- 0.2220 

0.1466 

0.2557 

- 0.3347 


0.0679 

- 

0.4227 

0.0012 

- 0.1450 

- 0.7998 

- 0.5402 

0.3123 

0.0532 

- 0.4019 

12 

0.0000 

- 

0.1350 

- 0.0016 

- 0.3108 

0.7750 

- 0.1569 

0.1307 

0.0841 

- 0.2253 


0.1181 

- 

0.1243 

- 0.0016 

0.6162 

- 1.5550 

- 0.1698 

0.1414 

0.0189 

- 0.1819 

13 

0.0000 

- 

0.1170 

- 0.0015 

0.0244 

0.0049 

- 0.1471 

0.1175 

0.0641 

- 0.1865 


0.1639 

- 

0.0639 

- 0.0009 

1.1130 

- 2.3520 

- 0.0906 

0.0737 

0.0071 

- 0.0914 

14 

0.0000 

- 

0.1046 

- 0.0016 

1.4280 

- 2.3960 

- 0.1408 

0.1040 

0.0581 

- 0.1579 


0.2246 

- 

0.0189 

- 0.0004- 

1.6570 

- 3.0610 

- 0.0301 

0.0240 

- 0.0057 

- 0.0239 

15 

0.0000 

- 

0.1323 

- 0.0066 

7.6720 

- 11.9200 

- 0.1912 

0.1158 

0.0816 

- 0.1780 


0.3068 


0.0349 

0.0019 

- 0.6730 

0.5485 

0.0449 

- 0.0246 

- 0.0363 

0.0485 

16 

0.0000 

- 

0.5012 

- 0.1355 

37.5100 

- 76.2100 

- 0.7038 

0.3409 

0.2381 

- 0.7026 


0.3412 


0.4944 

0.3420 

- 18.1200 

78.4600 

0.5483 

- 0.1452 

- 0.1522 

0.7367 

17 

0.0000 

- 

0.2754 

0.6667 

90.5100 

- 7.2760 

- 0.8031 

0,6909 

0.6666 

- 4.2170 


0.3434 


0.9189 

0.7194 

- 29,0500 

168.8000 

0.9950 

- 0.2849 

- 0.1838 

1.5530 

18 

0.0000 


0.7809 

0.0304 

- 89.6700 

120.3000 

1.1130 

- 0.5521 

- 0.7124 

1.0020 


0.3469 


1.3760 

- 0.8157 

-220.4000 

137.1000 

2.5630 

- 1.5400 

- 1.6360 

1.0260 

19 

0.0000 

- 

0,0097 

0.0065 

- 3.1110 

2.2820 

- 0.0285 

0.0242 

0.0258 

- 2.4050 


0.4312 


0.0115 

- 0.0080 

5.8500 

- 4.6220 

0,0254 

- 0.0282 

0.0224 

0.0314 

20 

0.0000 


0.3911 

- 0.1390 

63.6500 

- 22.3600 

1.1070 

- 0.3858 

- 1.0980 

0.3797 


0.4917 


0.2532 

- 0.0870 

10.1100 

- 4.3860 

0.7313 

- 0.2731 

- 0.6346 

0.2498 

21 

0.0000 


0.0027 

- 0.0026 

0.9412 

- 0.6377 

0.0001 

- 0.0066 

0.0044 

0.0151 


0,5654 


0.0027 

- 0.0019 

1.0500 

- 0.9433 

- 0.0011 

- 0.0299 

- 0.0375 

0.0538 

22 

0.0000 


0.2691 

- 0.0726 

18.2600 

- 5.4870 

1.2830 

- 0.3560 

- 1.4590 

0.4108 


0.6600 

- 

0.1539 

0.0343 

- 41.9300 

9.9780 

- 0.7248 

0.1372 

0.8950 

- 0.1667 

23 

0.0000 

- 

0.0263 

- 0.0005 

- 17.4500 

3.5250 

- 0.1362 

- 0.0131 

- 0.3378 

0.1425 


0.6677 


0.1667 

0.0404 

- 16.6400 

4.1580 

- 0.8415 

0.1809 

1,4090 

- 0.3044 1 

1 
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Table 6. Compensation Matrix [G(s)] in Cascaded 

Factor Form for the Space Shuttle Example. 




COMPENSATOR (1,1): GAIN 

= 1.0000 


COMPENSATOR COEFFICIENTS: 






ZA 

= 

.100000-01 

ZB 

= 

.744312 




ZC 

= 

1 .00000 

zo 

= 

1.02373 

ZE 

= 

8.16526 

ZC 

— 

1.00000 

ZD 

:= 

2.30558 

ZE 


2.01630 

ZC 

= 

1.00000 

ZD 

= 

1.22262 

ZE 

= 

4.44622 

ZC 


1.00000 

ZD 

= 

1.02373 

ZE 


8.16526 

PA 

= 

.133000-01 

PB 

= 

1.18973 




PC 

= 

1,00000 

PD 


4.36098 

PE 

= 

18.9002 

PC 


1.00000 

PD 

= 

4.36098 

PE 

~ 

18.9002 

PC 

= 

1 .00000 

PD 


2.38040 

PE 

= 

3,99765 

PC 

= 

1.00000 

PD 

= 

4.36098 

PE 

= 

18.9002 



COMPENSATOR (1,2): GAIN 

= .74500 


COMPENSATOR COEFFICIENTS 






ZA 

— 

1.00000 

ZB 

s 

-.996201 




ZA 


.100000-01 

ZB 

= 

1.01163 




ZC 

= 

1.00000 • 

ZD 


.582803 

ZE 

25 

8.16205 

ZC 

~ 

1.00000 

ZD 

= 

.656982 

ZE 

= 

10.4091 

ZC 

= 

1.00000 

ZD 

= 

.576605 

ZE 

■= 

2.04025 

ZC 


1.00000 

ZD 


1.20585 

ZE 

= 

4.44328 

PA 

= 

1.00000 

PB 

= 

• 50.0069 




PA 


.133300-01 

PB 

=: 

.983323 




PC 

= 

1.00000 

PD 

= 

5.20890 

PE 

= 

18.9033 

PC 


1.00000 

PD 


5.99161 

PE 

= 

25.0002 

PC 

= 

1.00000 

PD 

= 

5.32565 

PE 

= 

11.1114 

PC 


1.00000 

PD 

= 

2.39359 

PE 

= 

4.00060 



COMPENSATOR (1,3): GAIN 

= ,74500 


COMPENSATOR COEFFICIENTS 






ZA 


.000000 

ZB 

= 

51.0048 




ZA 

= 

.100000-01 

ZB 

= 

.869534 




ZC 


1.00000 

ZD 

= 

.583262 

ZE 

= 

8.16477 

ZC 

— 

1.00000 

ZD 


.657081 

ZE 


10.4118 

■ ZC 


1.00000 

ZO 

= 

.582806 

ZE 

= 

2.04267 

ZC 

= 

1.00000 

ZD 

=: 

1.21041 

ZE 


4.44575 

PA 

= 

1.00000 

PB 


50.0095 




PA 


.133000-01 

PB 

= 

1.10453 




PC 


1.00000 

PD 

= 

5.21469 

PE 

= 

18.9011 

PC 

= 

1.00000 

PD 

= 

5.99959' 

PE 

25 

24.9981 

PC 

=r 

1.00000 

PD 

i= 

5.33026 

PE 

= 

11.1092 

PC 

= 

1,00000 

PD 

= 

2.39180 

PE 

= 

3.99822 



COMPENSATOR (2,4): GAIN 

= 1.0000 


COMPENSATOR COEFFICIENTS 






ZA 

- 

.300000-01 

ZB 

= 

.995009 




ZC 


1.00000 

ZD 

:= 

.827221 

ZE 

= 

8.16655 

ZC 

= 

1.00000 

ZD 

= 

.992354 

ZE 

= 

6.25852 

ZC 

= 

1.00000 

ZD 

= 

2.00871 

ZE 

= 

1.56803 

ZC 

= 

1.00000 

ZD 

= 

1.20369 

ZE 

= 

4.45146 

PA 

= 

.400000-01 

PB 

= 

1.00117 




PC 

s 

1.00000 

•PD 

= 

4.34837 

PE 

= 

18.9029 

PC 

= 

1.00000 

P.D 

= 

4.34837 

PE 

= 

18.9029 

PC 

= 

1.00000 

PD 

= 

3.33705 

PE 

= 

11.1106 

PC 

= 

1,00000 

PD 


2.39557 

PE 


3.99576 
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Table 6. (Continued) Compensation Matrix [G(s)] in 
Cascaded Factor Form for the Space Shuttle. 




COMPENSATOR 

(1,1): GAIN = 

1.0000 


COMPENSATOR COEFFICIENTS 






ZA 


.100000-01 

ZB 

= 

1.00000 




ZC 

= 

1.00000 

ZD 

= 

1.00000 

ZE 

= 

8.1 6300 

ZC 

s 

1.00000 

ZD 

= 

2.28600 

ZE 

= 

2.01400 

ZC 


1.00000 

ZD 

= 

1.20000 

ZE 

= 

4.44400 

ZC 

= 

1.00000 

ZD 


1.00000 

ZE 

= 

8.16300 

PA 


.133000-01 

PB 

= 

1.00000 




PC 


1.00000 

PD 

= 

4.34800 

PE 

a 

18.9030 

PC 


1.00000 

PD 

= 

2.34800 

PE 

= 

18.9030 

PC 


1.00000 

PD 

= 

2.40000 

PE 

= 

4.00000 

PC 

= 

1.00000 

PD 

= 

4.34800 

PE 

a 

18.9030 



COMPENSATOR 

(1,2): GAIN = 

.74500 


COMPENSATOR COEFFICIENTS: 






ZA 

= 

1.00000 

ZB 

=S 

-1.00000 




ZA 

= 

.100000-01 

ZB 

= 

1 .00000 




ZC 


1.00000 

ZD 

= 

.571400 

ZE 

= 

8.16300 

ZC 

= 

1.00000 

ZD 

= 

.645200 

ZE 

a 

10.4100 

ZC 


1.00000 

ZD 


.571400 

ZE 

= 

2.04100 

ZC 

= 

1.00000 

ZD 

= 

1.20000 

ZE 


4.44400 

PA 

K 

.133000-01 

PB 

= 

50,0074 




•PA 

= 

1.00000 

PB 

= 

1.00000 




PC 


1.00000 

PD 

= 

5.21700 

PE 

= 

18.9030 

PC 

= 

1.00000 

PD 

= 

6.00000 

PE 

= 

25.0000 

PC 

= 

1 .00000 

PD 

= 

5.33300 

PE 

=r 

11.1111 

PC 

=: 

1 .00000 

PD 

= 

2.40000 

PE 

=r 

4.00000 



COMPENSATOR 

(1,3): GAIN = 

.74500 


COMPENSATOR COEFFICIENTS 






ZA 

= 

.000000 

ZB 

= 

51.0074 




ZA 

= 

.100000-01 

ZB 

= 

1.00000 




ZC 

=: 

1.00000 

ZD 


,571400 

ZE 

= 

8.16300 

ZC 


1.00000 

ZD 

= 

.645200 

ZE 

a 

10.4100 

ZC 

= 

1.00000 

ZD 

= 

.571400 

ZE 

= 

2.04100 

ZC 


1.00000 

•ZD 

r: 

1.20000 

ZE 

= 

4.44400 

PA 

= 

1.00000 

PB 

= 

5010074 




PA 


.133000-01 

PB 


1.00000 




PC 

s; 

1.00000 

PD 

= 

5.21700 

PE 

a 

18.9030 

PC 

s 

1.00000 

PD 

= 

6.00000 

PE 

a 

25.0000 

PC 

= 

1.00000 

PD 


5.33300 

PE 

a 

11.1110 

PC 

= 

1.00000 

PD 

= 

2.40000 

PE 

= 

4.00000 



COMPENSATOR 

(2,4): GAIN = 

.74500 


COMPENSATOR COEFFICIENTS 






ZA 


.300000-01 

ZB 

= 

1.00000 




ZC 


1.00000 

ZD 

= 

.857000 

ZE 

a 

8.16300 

ZC 

=r 

1.00000 

ZD 

= 

1.00000 

ZE 

= 

6.25000 

ZC 


1.00000 

ZD 

a 

2.00000 

ZE 

= 

1.56300 

ZC 

= 

1.00000 

ZD 

= 

1.20000 

ZE 

a 

4.44400' 

PA 

r= 

.400000-01 

PB 

= 

1,00000 




PC 

= 

1.00000 

PD 


4.34800 

PE 

a 

18,9030 

PC 

= 

1,00000 

PD 

= 

4.34800 

PE 

a 

18.9030 

PC 

r 

1.00000 

PD 

- 

3.33300 

PE 

a 

11.1110 

PC 


1.00000 

PD 


2.40000 

PE 


4,00000 

COMPENSATORS HAVING ZERO 

CONTRIBUTION: (1.4) 

(2,1); 

(2.2); (2,3) 
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Table 7. System Specifications for the Space Shuttle Example. 


Margin 

Number 

Margin 

Value 

Complex Frequency 
RE[s] IM[s] 

Desired 

Margin 

Margin 

Type 

Active 

List 


A. 

Subsystem No. 1, 

Iteration 

No. 0 

• 

1 

0.4163 

0.000 

0.0896 

0.60 

GM 

Yes 

2 

53.7.400 

0.000 

0.0240 

30.00 

PM 

No, 

.3 

i27.5000 

0.000 

0.0603 

■'30.00 

PM 

No 

4 

24.2100 

0.000 

0.0679 

•30.00 

PM 

Yes 

5 

0.0074 ■ 

0.000 

0.3412 

0.10 

AM 

No 

6 

0.0289 

0.000 

0.3469 

0.10 

AM 

No 

7 

8 

0.0097 

0.000 

0.6600 

0.10 

AM 

No 


Subsystem No. 2, 

Iteration 

No. 0 


1 

0.5941 

0 . 000 

0.0896 

.0.60 

GM 

Yes 

2 

36.5000 

0.000 

0.0321 

30.00 

PM 

No 

3 

0.0885 

0.000 

0.3434 

0.10 

AM 

Yes 

4 

0.0076 

0.000 

0.6600 

0.10 

AM 

No 


B. 

The Improved System Specifications 
Subsystem No. 1,. Iteration No. 5 


1 

. 0.6120 

0.000 

• 0.0927 

0.60 

GM 

No 

2 

38.1800 

0.000 

0.0222 

30.00 

PM 

No 

3 

143.9000 

0.000 

0.0616 

30.00 

PM 

No 

4 

32.0100 

0.000 

0.0674 

30.00 

PM- 

No 

5 

0.0307 

0.000 

0.3469 

0.10 

AM 

No 

6 

0.0098 

0.000 

0.4917 

0.10 

AM 

No 

7 

0.0041 

0.000 

0.6600 

0.10 

AM 

No 



Subsystem No, 2, 

Iteration 

No. .5 


1 

0.6075 

0.000 

0..0896 

0.60 

GM 

No 

2 

31.3600 

0.000 

0.0339 

30.00 

PM 

No 

3 

0.0845 

0.000 

-0.3434 

0.10 

AM 

No 

4 

0.0076 

0.000 

0.6600 

0.10 

-AM 

No 




CHAPTER V 


CONCLUSIONS AND RECOMMENDATIONS 
Conclusions 

Because of the complexity of technology and control laws, the 
design of modern control systems has become Increasingly compli- 
cated. In this exposition, the theory and associated numerical 
technique for achieving a computer-aided compensation design improve- 
ment algorithm for the multivariable control system have been pre- 
sented. The technique developed is applicable to linear, time- 
invariant systems possessing multiple input, multiple output status 
whose plant characteristics are described by discrete open-loop 
frequency response data. The compensation matrix is entered as 
transfer functions of cascaded first and second order polynomials. 

The method was designed using a strict constraint algorithm to 
alleviate the inherent problems generally associated with soft con- 
straint cost functionals. The objective of the Compensator Improve- 
ment Program is to modify in an iterative manner the free parameters 
of the compensation yielding a system that satisfies specified 
frequency response properties. The computer coding in the Fortran 
IV language has been included and the practical utility of the 
program illustrated with space related examples. 

Chapter I contains a literature survey of the previous 
research on the automatic design problem based on theory developed 



60 


by classical design methods. Each of the aforementioned techniques 
has been ascertained successful for the author's specified control 
area, generally restricted to a single control input system; 
however, the Compensator Improvement Algorithm [7] appeared most 
readily applicable and available for extension to the multivariable 
control case with the objective of obtaining a suboptimal solution 
to the specified constraints of the control design. 

The theoretical concepts of the design algorithm are developed 
in Chapter II. The mathematical derivations associated with the 
algorithm in determining the necessary closed-loop frequency 
response, gradient vectors, and hence, the directional vectors , 
have been deduced by exact means. 

In Chapter III the objective of the schematic algorithm and the 
computational flow diagram were introduced. This chapter contains 
an explanation of the flow diagram referring to the theory and 
synopsis of the subprograms in Appendix A. 

Pragmatic examples illustrate the effectiveness of the 
algorithm in Chapter IV'. Three space related examples were pre- 
sented: Cl) a dual input, dual output system exemplifying no 

coupling; (2) a dual input, dual output system with coupling; 

C3) a Space Shuttle example with dual control inputs and four out- 
puts. The results herein verify the utility of the Compensator 
Improvement Program as a practical method of compensation improve- 


ment. 
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In conclusion, this exposition has demonstrated that the classi- 
cal control theory is amenable to systems with multivariable charac- 
teristics. An important benefit derived from the use of the frequency 
domain for linear time-invariant systems is the intuition it provides 
in determining the soundness of a system. The digital computer has 
made possible the application of classical techniques to the optimal 
design problem. The ultimate contribution of this research effort 
is the development and implementation of an algorithm to enhance 
compensator design of systems possessing multivariable control char- 
acteristics. 


Recommendations 

The employment of any digital computer algorithm as an aid in 
the aggregate design process is perhaps ^s much an art as- the design 
process itself.' The use of the computer-aided design program may 
free the engineer from many burdensome and time-consuming calculations, 
but it is the engineer who in essence must .provide the framework in 
which to enter the compensation in order to achieve the desired control 
law. This then is perhaps the greatest limitation of any computer- 
aided control design; that is, there is no supplanting the awareness 
and judgement of an experienced- control engineer. 

More realistically, however, the Compensator Improvement Algorithm 
does possess minor limitations which could be reconciled. In particu- 
lar, 

1. CIP should be given the option of accepting either 
discrete frequency response data or transfer function 
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information. With the transfer function option, CIP 
would calculate its own frequency response data. 

2. Modifications should be made for CIP to accept 
prefilters; these filters would be user designated 
and not altered by the program. 

3. CIP should be given the option of accepting sampled- 
data systems without the necessity of the engineer 
converting compensation into the W-plane; possibly 
CIP could be further extended to accept multirate 
sampling problems, 

4. The practical utility of CIP could be extended by 

rendering the algorithm capable of producing a two 
phase optimization program: in particular, the 

present version of the algorithm would yield a de- 
sign to meet a set of design objectives producing a 
feasible solution while continuing to satisfy the 
desired specifications. For the second phase, per- 
haps a gradient projection technique could be uti- 
lized in optimizing the necessary cost function. 

In essence the major limitation of the CIP algorithm is its 
restriction to linear, time-invariant control systems. The aspect: 
of extending the work to nonlinear systems have not been contempla- 
ted; this omission is regrettable since the occurrence of system 
uncertainty is always a possibility. In regard to the restriction 
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of time- invariance, the present CIP techniques are hot amenable 
to time-varying systems. 

Mississippi State University 

Mississippi State, Mississippi 39762 
July 25, 1977 



APPENDICES 



APPENDIX A 


SUBPROGRAM SUMMARIES OF THE COMPENSATOR 
IMPROVEMENT PROGRAM 

Introduction 

It is the objective of this Appendix to provide the basic con- 
cepts in theory and/or programming techniques incorporated within 
each subroutine. With the enclosed information any efforts made in 
adaptations or modifications for solving related problems should be 
reduced significantly. The subprograms are presented in alphabetical 
order for easy reference. 
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Subroutine - ADOPTS 

In an effort to minimize the input data storage required and 
the corresponding computer time used in manipulating extensive data, 
the subprogram ADDPTSI9] generates additional frequency data. In 
particular, if the spacing of the original response data in the 
neighborhood of a critical point in the relative stability region 
becomes too large, this subprogram interpolates the given data in 
this neighborhood yielding a more accurate stability margin. This 
design philosophy is based on the continuous nature of the frequency 
response over the complete range of frequencies. The added frequency 
points are obtained in accordance with the routine INTER, an inter- 
polation algorithm; a log type of interpolation is used in 
determining all magnitudes; whereas, phases are calculated by linear 
interpolation. 

This subroutine also requires the .routines CRT and EVAL for 
updating the frequency response at the data points. The routine 
DELETE is used in conjunction with ADDPTS to retain only the origi- 
nal data at each new iteration. 

The following variables are designated for this subprogram: 

Subprogram Variables 

■Input Variables ; 

KPOINT - An integer variable used to denote the current number 
of data points. 

KIN - An integer variable that denotes the number of inputs 


to the controller. 
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KOUT 

NB 

NM(I) 

STBM(I) 

KPTS(I,J) 

CT(I,J) 

G(I,J,K) 

GC(I,J,K) 

T(I,J,K) 

OMEGA(I) 

KPTMAX 


- Integer variable denoting the number of controller 
outputs . 

- An integer used as a counter representing the starting 
number of the margins to be investigated. 

- An integer array representing the number of margins 
investigated for each subsystem. 

- A one dimensional real array containing the values of 
the stability margins. 

- A two dimensional integer array containing the fre- 
quency numbers of the margins in accordance with a 
particular subsystem. 

- A two dimensional complex array containing the overall 
frequency response of the system. 

- A three dimensional complex array containing the origi- 
nal discrete data frequency response of the plant 
system. 

~ A complex three dimensional array storing the compen- 
sator response evaluated at the specified frequency 
points. . 

- A complex three dimensional array in which the transfer 
response [GC]*[G] is stored. 

- A one dimensional complex^ array containing, the discrete 
frequency points. 

- An integer denoting the maximtim number of discrete 
frequency points allowable on a single iteration. 
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KGOBAK - An integer variable used to denote whether frequency 

points were added. 

NIT (I) - An integer array denoting the number of Inactive 

margins for each subsystem. 

KINACT(I,J) - A two dimensional array containing the integer numbers 
corresponding to the frequencies of the inactive or 
satisfied margins for each subsystem. 

NML(I) - An integer array denoting the number of active margins 

detected per subsystem. 

KACT(I,J) - A two dimensional integer array denoting the frequency 
data numbers of the active or unsatisfied margins for 
each subsystem. 

KPOLD - An integer denoting the number of data points on the 

last iteration. 

KOLD(I) - A one dimensional integer array containing the previous 

data points. 

KSYM - An integer variable used to reference the frequency 

response of the particular subsystem being manipulated. 

The following transient variables are used in the auxiliary 

subroutines CRT and EVAL and are not affected directly by this sub- 
program; for more information regarding these variables refer to 

the respective subroutine synopsis. 

CRT: C1,CI,W0RKI 


EVAL; ZA,ZB,ZC,ZD,ZE,PA,PB,PC,PD,PE, 
N1,N2,M1,M2,GAIN,K0NT,A,B,C»D,E 
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Output Variables ; 

The following output variables are defined as their respective 
input counterparts, but as outputs have been updated to include the 
newly generated data points: G(I,J,K), GC(I,J,K), T(I,J,K), KPOINT, 
OMEGA(I), KINACT(I,J), KACT(I,J), KOLD(I). 
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Subroutine CHANGE 

The subprogram CHANGE is designed to change the individual com- 
pensator coefficients in an orderly manner to induce an improved 
system. The change in the compensator coefficients is made in 
accordance with the directional vector of the particular iteration. 

Note that the compensator elements are described by the transfer 
functions in the form of cascaded first and second order factors, 
that is, 

N1 N2 

n (ZA. + ZB.s) H (ZC, + ZD.s + ZE.S ) 

G(s) = (GAIN) ^ — • (A-1) 

Ml M2 

n (PA. + PB.s) n (PC. + PD.s + PE.s^) 

i=l ^ ^ j=l J J 

Recall that the ultimate goal of the Compensator Improvement 
Program is the design of. compensation so that the measurements of 
the system performance are equal to or better than the system speci- 
fications. The design of the compensators can be expressed mathe- 
matically as the strict constraint problem: 

T 

Determine x 

(A-2) 

subject to g^(x^) ^ b^, i = l,...,m 

where x is a vector of the n compensator coefficients. The functions 
g^(x^) contain measurements of the system performance in terms of 
the compensator coefficients; thus, these functions represent the 
stability and attenuation margins, as well as any constraints on the 
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compensator coefficients. The constants are the system specifi- 
cations . 

Assuming that for equation (A-2) the trial solution vector at 
T 

the kth iteration is , then a trial solution vector of a possible 
improved solution at the (k + l)th iteration is 

+ h[VG] a (A-3) 

where [VG] is the Jacobian matrix evaluated at and consists of 
all the active constraints, i.e., the functions g.(xv) < b, . The 

1 K X 

scalar h is a normalized step size. • The vector ^ is calculated as 

a = [VG^VG]“^£ . (A-4) 

where ^ is a vector of weighting constants initially set at unity. 
The routine DIRVEC determines the directional vector ^ , where 

d = [VG]a (A-5) 

Thus the subprogram CHANGE ultimately applies the elements of this 
directional vector to its corresponding compensator coefficients 
weighted by the step size h. 

' The program variables- are defined' as follows: 

Subprogram Variables 

Input Variables ; 

ZA(I,J,K) - A real three dimensional array representing the con- 
stant terms of the first order factors of the 
compensation polynomials. 
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ZB(I,J,K) - A real three dimensional array containing the first 
order factor coefficients of s. 

ZC(I,J,K) - A real three dimensional array representing the con- 
stant terms of the second order factors. 

ZD(I,J,K) - A three dimensional real array containing the s coef- 
ficients of the second order factors. 

ZE(I,J,K) - A real three dimensional array containing the s 
coefficients of the second order factors. 

N1(I,J) - An integer array that denotes the number of cascaded 

first order factors. 

N2(I,J) - An integer array that denotes the number of cascaded 

second order factors. ■ 

DV(I) - A real one dimensional array representing the 

directional vector ^ in equation (A-5). 

DEL - A real variable that- denotes the step size h in (A-3) . 

KKK - An integer used to count the nximber of elements in 

the directional vector. • 

KIN - An integer denoting the number of controller inputs 

or sensed states. 

KOUT - An integer denoting the number of controller outputs. 

A,B,C - Parameter variables used to dimension the arrays by 

the number of maximum allowable cascaded factors. 

Output Variables : 

The following output variables are defined in the same manner as 

their respective input variables, but have been updated incremental- 
ly in accordance with the directional vector and step size: ZA(I, J,K) , 

ZB(I,J,K), ZC(I,J,K), ZD(I,J,K), ZE(I,J,K). 
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Subroutine CRT 

The subprogram CRT has two major functions as denoted by the 
input variable KEY. If KEY has been set to unity, the program 
calculates the closed-loop frequency response [C(s)] in terms of 
the input vector, that is, 

[C(s)j = [G(s)][P(s)]{l + [H(s)][G(s)][P(s)]}"^[R(s)] . (A-6) 

Thus the open-loop frequency response of the kth system can be 
obtained by setting the Rth diagonal element of [H(s)] to zero, and 
by setting all the elements of [R(s)] to zero except the kth element 
which' is set to unity; the result is the frequency response between 
the kth input and the outputs when the kth loop is open. 

If KEY is entered as zero, the program is designed to aid the 
subprogram PARTAL in. determining the partial term 3Cj^(s) /9G^^. . 
Actually, this term can be evaluated by the equation 

■41^ = -IG(s)][P(s)3{I + rH<s)][G(s)HP(s)3r' [H(s)3 • 

ij ij 

[P(s)3{i + EH(s)3[g(s)3[p(s)3}“'[RCs)3 

+ [P(s)3{i + [h(s)3EG(s)3[p(s)3}"'[R(s)3 (a-7) 

ij 

which bears semblance to equation (A-6). Then, 3C, (s)/3G.. is the 

K 13 

kth element of (A-7) with the conditions aforementioned. 

Other routines used by this subprogram are the complex matrix 
inversion and multiplication programs, MATING and MATMUL 
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respectively. The variables are defined in the following list: 

CRT Subprogram Variables 

Input Variables ; 

KEY - A switch variable which designates the subprogram 

mode: if KEY is set to unity, the response of (A-6) is 
returned; if KEY is set to zero, the partial term of 
(A-7) is calculated and returned. 

K - An integer variable denoting specific frequency point 

under consideration. 

KSTM - An integer variable representing the system under con- 

sideration, that is, k in the previous discussion. 

T(I,J) ■ - A two dimensional complex array of the transfer response 

[G(s)] • [P(s)]. 

KIN - An integer denoting the number of control inputs. 

KODT - An Integer denoting the number of controller outputs. 

P(I,J) - A two dimensional complex array of the plant response. 

II - An integer denoting the control input index of G.. in 

(A-7), 

J1 - An integer denoting the control output index of the 

term in equation (A-7) . 

Cl(l,J) - A two dimensional internal complex array containing the 
total frequency response at a particular frequency 
point . 

CI(I,J) - An internal two dimensional array storing the inverse 
of the complex frequency response Cl(l,J), 
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W0RKI(I,J) - A two dimensional complex array for internal subroutine 
use. 

A,B,D - Variables denoting dimensxon axxocations. 

♦ * 

Output Variables ; 

CT(I) - A one dimensional complex array representing the 

response [C(s)] in equation (A-6). 

PCG - A complex variable representing the partial term 

8C, (s)/3G.. in the previous discussion. 

1C xi 
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Subroutine DELETE 

The subprogram DELETE is designed in conjunction with the ADDPTS 
routine in an effort to save computer storage and time in manipu- 
lating unnecessary data. In particular, this subprogram deletes the 
extra frequency points and their corresponding response terms gene- 
rated by the ADDPTS routine when the maximum number of allowable 
points has been generated; for accuracy, the original input data is 
always retained. 

The" routine uses the following variables; 


Subprogram Variables 


Input Variables ; 


KPOINT 

KIN 

KOUT 

OMEGA (I) 
G(I,J,K) 

NIT(I) 

KINACT(I,J) 


An integer counter to denote the number of frequency 
points . 

An Integer denoting the number of controller inputs 

V 

or the measured states. 

An integer variable representing the number of con- 
troller outputs or plant inputs from the controller. 
A one dimensional complex array of. frequency terms. 

A three dimensional complex array consisting of the 
frequency response describing the plant system. 

An integer array denoting the number of inactive or 
satisfied margins detected. 

A two dimensional integer array of the frequency 
numbers of the corresponding inactive margins for 
each subsystem. 
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NML(l) 

KACT(I,J) 

ITER 

KPOLB 

KOLD(I) 

KNEW(I) 

A,B,C,D,E 


- An integer array denoting the total number of active 
or unsatisfied margins detected per subsystem. 

- A two dimensional integer array of the frequency 
data of the active margins. 

- An integer representing the present iteration. 

- An integer denoting the total number of frequency 
points. 

- An integer array of the original frequency points 
reference numbers. 

- An integer array of the generated frequency numbers. 

- Dimension allocations; set by parameter statement 
in the main program. 


Output Variables ; 

The following output variables correspond to their respective 

I 

input variables, but have been updated by the subprogram deleting 
the extra data; OMEGA (I) , G(I,J,K), KINACT(I,J), KACT(I,J), KOLD(I). 
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Subroutine DIRVEC 

The objective of the subprogram DIRVEC is to calculate the 
directional vector of the constraint improvement algorithm. The 
directional vector _d is calculated as 

d = [VG] a (A-8) 

where VG is a matrix of order (n,m) whose columns consist of the 
gradients of the active constraints. The m component column vector 
is determined by 

a = [VgV]“^ ^ (A-9) 

where c is a m component column vector whose elements are all 
positive. The order m corresponds to the number of compensator co- 
efficients. Auxiliary subprograms include the matrix inversion for 
real variables MATINV and the matrix multiplication MATMUL routines. 
Definitions of the input and output variables are as follows: 

Subprogram Variables 

Input Variables ; 

G(I,J) - A two dimensional array whose coltunns are comprised 

of the gradients of the active constraints, 

NM - Integer value used to designate the number of col^lmns 

in G, that is, the number of active constraints. 

KPARC - Integer variable that represents the number of rows 

in G, that is, the number of compensator coefficients. 
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WEIGHT (I) - A real one dimensional array that contains the 

column matrix jc of equation (A-9). 

E,F,Z,H - Parameter variables for allocating dimension storage. 

The transient variables, A1,AI, and WORKR, are used in the 
auxiliary subprograms MATIKV and MATMUL and are not affected direct- 
ly in this routine. 

Output Variables ! 

DV (I) - A real one dimensional array which corresponds to the 

directional vector d in (A— 8) . 



80 


Subroutine EVAL 

The subprogram EVAL is designed to evaluate the compensator 
matrix [G(s)] at each discrete frequency point. Recall the compen- 
sation elements are polynomials of the form of cascaded first and 
second order factors, that is, 


G(s) = (GAIN) 


II (ZA. + ZB.s) n (ZC. + ZD.s + ZE.s ) 
!=i . ^ i 3 3 3 ' 


n (PA. + PB.s) n (PC, + PD.s + PE.s^) 
t-1 ^ ^ 3=1 J ^ ^ 


The subprogram utilizes the polynomial evaluation routine POLEV. 

This subroutine also calculates the transfer relation [G(s)]* 
[P(s)l for each frequency point using the matrix multiplication 
program MATMUL . 

The subprogram uses the following variables: 


Subprogram Variables 

Input Variables ; 

XF -A complex variable denoting the discrete frequency 

point . 

G(I,J,K) - A three dimensional complex array containing the 

original discrete data frequency response of the plant. 
K - An integer variable used to denote the current 

number of the frequency data point. 

KIN - An integer variable denoting the number of Inputs 


to the controller. 
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KOUT 

ZA(I,J,K) 

ZB(I,J,K) 

ZC(I,J,K) 

ZD(I,J,K) 

ZE(I,J,K) 

PA(I,J,K) 

PB(I,J,K) 

PC(I,J,K) 

PD(I,J,K) 

PE(I,J,K) 

N1(I,J) 

N2(I,J) 


- An integer that denotes the number of controller out- 
puts, 

- A real three dimensional array containing the constant 
terms of the first order factors of the numerator. 

. - A real three dimensional array containing the first 
order factor coefficients of s in the numerator. 

- A real three dimensional array representing the con- 
stant terms of the second order numerator factors. 

- A real three dimensional array , containing the s coef- 
ficients of the second order factors of the numerator. 

- A real three dimensional array storing the s coeffi- 
cients of the second order factors -in the numerator, 

- A real three dimensional array containing the constant 
terms of the first order denominator factors. 

- A .real three dimensional array containing the first 
order denominator factor coefficients of s, 

- A real three dimensional array representing the con- 
stant terms of the second order denominator factors, 

- A real three dimensional array containing the s coeffi- 
cients of the second order factors in the denominator. 

- A real three dimensional array storing the s^ coeffi- 
cients of the second order denominator factors, 

- An Integer array that denotes the number of cascaded 
first order factors in the numerator. 

- An integer array that denotes the number of second 
order cascaded factors in the numerator. 
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- An integer array denoting the number of cascaded first 
order factors in the denominator. 

M2(I,J) - An integer array that denotes the number of second 

order cascaded factors in the denominator. 

GAIN (I, J) - A two dimensional real array containing the DC gain 

of each compensator polynomial. 

KONT(I,J) - A two dimensional integer array designating whether 
the DC gain for the particular channel. is allowed to 
vary: if KONT is unity, the gain may vary; if KONT 

is two, the gain is- not allowed to vary. 

A,B,C,D - These -variables are defined by a parameter statement 

in the main program and designate maximum storage al- 
locations in regards to the order of the arrays. 

Output Variables ; 

GC(I,J,K) - A three dimensional complex array storing the compen- 
sator response evaluated at specified frequency points. 

T(I,J,K) - A complex three dimensional array in which the transfer 
response [GC] * [G] is stored. 
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Subroutine GAINMG 

The purpose of the subprogram GAINMG is to locate and calculate 
the gain margins of a system represented by a discrete open-loop 
frequency response. With f^ as the ith frequency specified, the 
corresponding complex frequency response can be represented in real 
and imaginary terms as GR^ and GI^ respectively. Thus the following 
sequence can be formed to detect the occurrence of a gain margin; 

U. = GI. • GI. . . (A-11) 

1 1 1-1 

A gain margin is located whenever becomes either negative or 

zero. The frequency number of the gain margin is taken as i or i-1 
depending on whether [Gl^j > or _< > the gain 

margin is calculated as 

STBM = ll. + GR^ + jGIj^l , (A-12) 

where k is either i or i-1. 

The following variables are designated for the subprogram; 
Subprogram Variables 

Input Variables ; 

OMEGA(I) - A complex one dimensional array that contains the 

specified frequencies in ascending order for 
describing the system. 

GTOTAL(I) - A complex one dimensional array containing the com- 
pensated open-loop frequency response for the Ith 
frequency. 
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KPOINT - An integer number of frequency points used to des- 

cribe the open-loop frequency response of the system. 

FQMIN - The lowest frequency for which gain margins are 

detected. 

FQMAX - The largest frequency for which the gain margins are 

to be detected. 

NM - The integer used as a counter for the number of 

margins located. For example, assume HM is initially 
2 and this program locates 3 margins; these margins 
would be labelled as margins 3, 4, and 5' respectively. 

Output Variables : 

NM - This is the number that designates the last gain 

margin found. 

KPTS(I) - A one dimensional integer array that contains the 

frequency members of the margins found. 

STBM(I) - A one dimensional real array that contains the margin 

values corresponding to the frequency pointer KPTS, 
These margins are measured in terms of distances 
from the point (-1 + jO) in the complex GH(jw) 


plane . 
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Subroutine IISITER 

INTER is a subprogram designed in conjunction with the 
ADOPTS routine to interpolate specified input data, thereby 
generating new data. The algorithm is designed to yield a log type 
of interpolation in determining all magnitudes; whereas, phases are 
calculated by a linear interpolation scheme. 

The variables are defined in the listing: 

Subprogram Variables 

Input Variables ; 

S - A complex number consisting of the lower bound of the 

quantity to be interpolated. 

T - A complex number consisting of the upper quantity bound 

of the interpolation. 

Output Variable : 

R - A complex number representing the resultant of the 


interpolation. 
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Subroutine MATING 

Subroutine MATING determines the inverse of a matrix of complex 
elements by the Gauss-Jordon reduction method. It is assumed that 
no diagonal elements of the original matrix are zero. If in applying 
the reduction procedure the magnitude of the ith element of the ith 
pivot row is of magnitude less than 1.0 x 10 , the inverse matrix 

is assumed nonexistent. 

The input » output variables are defined as: 

Subprogram Variables 

Input Variables i 

XX(I,J) ■ - A complex two dimensional square array whose inverse 
is desired. 

N - An integer denoting the number of rows and columns in 

matrix XX(I,J). 

X(I,J) - A complex array containing the generated augmented 
matrix. 

A,B,C - Parameter variables denoting storage allocation. 

Output Variables : 

YY(I,J) - A complex two dimensional array that contains the 
inverse of the XX(I,J) array. 

lER - The error code of the subprogram. If lER is zero, 

no error was incurred; if lER is unity, the matrix 
is assumed to be singular. 
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Subroutine MATINV 

The subprogram MATINV uses the Gauss-Jordon reduction method 
in determining the inverse of a matrix of real elements. The 
procedure and program variables of this routine are defined in the 
same manner as in the subprogram MATING but with application to the 
real matrix XX (I, J) and its inverse YY(I,J) composed of real 
elements only. 
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This subprogram is designed to operate on either real or com- 
plex matrices as specified by the input variable NC. The input, 

output variables are designated as follows: 

Subprogram Variables 

Input Variables : 

AC(I,J) - A complex two dimensional array representing the 

matrix A as aforementioned when the subprogram is used 
in the complex mode. 

BC(I,J) - A complex two dimensional array representing the 
aforementioned matrix B. 

AR(I,J) - A real two dimensional array corresponding to matrix A 
in the previous discussion when the subprogram is 
used in real term mode only. 

BR(I,J) - A real two dimensional array corresponding to a real 
matrix B. 

N - An integer variable denoting the number of rows in the 

matrix A. 

L - An integer variable denoting the number of columns in A. 
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111 - An integer variable denoting the number of rows in 

matrix B in the above discussion. 

M - An integer variable that denotes the number of columns 

in matrix B. 

ND - An integer variable used to designate the proper 

storage placement when multiplying three dimensional 
matrices; for the two dimensional case set ND to unity. 

NC - An integer which designates whether the program is to 

multiply real or complex matrices. If NC is zero, the' 
complex matrices AC, BC, CC are used; if NC is unity, 
the real matrices AR, BR, CR are manipulated. 

A,B,C, - These variables are defined by a parameter statement in 
the main program and designate maximum storage capabili- 
ties in regard to the order of the matrices. 

Output Variables ; 

CC(I,J,ND) - A complex three dimensional array that contains the 
complex elements c^^ in equation (A-13) . 

CR(I,J,ND) - A real three dimensional array containing the resultant 
matrix C when the subprogram is in real mode. 
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• Subroutine NYQUIST 

The subprogram NYQUIST determines the number of closed-loop 
poles of a closed-loop system inside a certain enclosed contour of 
the s-plane. The number of closed-loop poles within the contour is 
from the relation 

Z = P + N , (A-14) 

where Z is the number of closed-loop poles, i.e., the number of 
roots of the characteristic equation inside the contour; p is the 
number of open- loop poles within the contour, and N is the algebraic 
sum of the encirclements around the (-1 +10) point by the frequency 
response. Note the encirclements around the (-1 + jO) point are 
assvimed positive if clockwise and negative if counterclockwise. 

The NYQUIST encirclements are counted by application of the 
equation: 

N = INTEGER ^ + I {ANGLE[1 + G(s.)] 

( ^ ^ i=2 ^ 

- ANGLE tl + G(s^_^)]}/7r| CA-15) 

where the operator INTEGER yields only the integer portion of the 
calculation in brackets, and' P is the number of poles on the con- 
tour. It is assumed that the contour is symmetric with respect to 
the real axis in the s-plane; hence the frequency response G(.s) is 
symmetric about the real axis of the G(s) -plane. Thus in 
equation (A-15), only frequency points on that portion of the 
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contour in the upper half-plane are used in the evaluation. The 
fraction i 1/8 is used to account for the fact that the frequencies 
and Sjj , the first and last frequency points respectively, may not 
actually he on the real axis. It is assumed that the points are 
chosen so that the sum of the angles is within tt/ 4 radians of the 
correct value. The positive and negative signs are chosen in agree- 
ment with the the sign of the summation respectively. 

The input, outpub variables are defined as follows: 

Subprogram Variables 

Input Variables : 

N - An integer variable that denotes the number of frequency 

response points supplied by the user, i.e,, M in (A-15). 
G(I) - A complex one dimensional array containing the frequency 

response G(s) where s is chosen along the contour. 

NRHP - An integer number of the open-loop poles inside the 

contour, i.e., P in the previous discussion. 

NCON - An integer number of open-loop poles located on the 

contour . 

FMAX - A real variable denoting the maximum frequency range. 

P(I) - A complex one dimensional, array containing the 

specified frequency points. 

Output Variables: 

NCIRL - An integer representing the number of encirclements 

around the (-1 + jO) point. 

NZ - An integer of the number of poles of the closed-loop 

system inside the contour. 
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Subroutine OUTPT 

The purpose of this subprogram is to output certain information 
at various stages of the main program. Three areas of information 
available for output are: Compensator Information, Frequency 

Response Information, and Stability Margin Data. 

If the user desires data on the compensation matrix, the. 
selector variable N is set to zero and the program outputs the 
compensator values at the last iteration. 

By setting the selector variable to unity, the overall frequency 
response data is the output. 

For investigating the stability margins, the selector variable 
is set to two. The corresponding output data includes: 

1. The margin numbers 

2. The frequency where each margin occurs 

3. The value of each margin 

4. The desired value of each margin 

5. The type of margin, i.e., phase margin (P) , gain margin 
(G) , stability margin (S) , or attenuation margin (A) 

The directional vector at the last iteration. 


6 . 
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Subroutine PARTAL 

The Compensator Improvement Program is given the ability to 
determine the necessary partial derivatives of the various active 
margins, assuming a frequency response along a general contour, by 
the subprogram PARTAL. In order to develop this subroutine, the 

i 

theoretical derivations of Section 3 of Chapter II, were implemented 
to yield these partial derivatives. 

Recall that in designing compensation, CIP searches the compen- 
sated frequency response C(jw) over specified ranges of frequency to 
determine which margins do not satisfy the desired values. For the 
unsatisfied or active margins, the design specifications are con- 
verted to distances between certain critical points of the open-loop 
frequency response and particular points in the corresponding complex 
plane; that is, the typical objective function for the kth open- loop 
system is 

d = {[A + Cj^(jw)][A + Cj^(ju)]*}^ (A-16) 

where Cj^(jto) is the kth system's compensated response and A is the 
point in the complex plane from which the specification is measured. 
Generally, point A is chosen as the (-1.0 + jO.O) point for stability 
margins and (0.0 + jO.O) for attenuation margins. 

PARTAL determines the gradients of these design objectives with 
respect to the compensator coefficients of the controller. Thus, the 
partial derivatives of d with respect to some parameter w is 
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9C, (jw) ) 

[A + C^(jo))]* ^3^— > - d . (A-17) 

Evaluation of (A-17) depends largely on determining 9C^(jw)/3w 
accurately. Using tlie chain rule- this becomes 



3C^(jw) 

3w 


3C, (joJ) 3G. . 

XI 

3G.. ‘ 3w 


(A-18) 


where G is the element of the controller [G(s)] in which the free 

ij 

parameter w appears. 

As derived in Chapter II, the first term in equation (A-18) 
may be evaluated by taking the kth element of 

= -[G(s)]-[P(s)]{l + [G(s)][P(s)][H(s)]}"^[H(s)] ^ 

ii 


[P(s)]{l+ tH(s)][G(s)][P(s)ll ^[R(s)l 

(A-19) 

+ [P(s)]{l + [H(s)][G(s)]tP(s)]r'[R(s)] 


where the kth diagonal element of [H(s)] is set to zero and all the 

elements of [R(s)] are set to zero except the kth which is set to 

unity. Note also that 3 [G(s)]/ 3G. . is a zero matrix except for 

the (ij)th element which is unity. 

The second term in (A-18), 3G../3w , is derived in Chapter II. 

^3 

Assuming the (ij)th element of [G(s)] is composed of a cascaded 
arrangement of transfer functions, i.e.. 
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G . , (s) 


K 

n 

k=l 



(A- 20) 


where K is the number of cascaded elements. The X-th cascaded 
element of the (ij)th compensator of [G(s)] has the general form 


G. . 


(s) 


M 

m=0 

L 


n=0 


(A-21) 


Then, the free parameters of this element are the x’s and y’s. If 
w in (A-17) is the pth numerator coefficient of (A-21) , then 


3G. Cs) 


+s^ 


3x 


ijS-p 




V m 

Z X. . „ s 
^ ij£m 
m=0 / 


(A-22) 


Similarly, if w represents the pth denominator coefficient of (A-21) , 
then 


9G. .(s) 
^^ij£p 


-s" 


G..(s) 


N 


^ ^iUn" 


,n=0 


(A-23) 


Equations (A-16) , (A-18) , (A-22) , and (A-23) indicate how the needed 
partial derivatives can be calculated. The subprogram PARTAL imple- 
ments these equations including the necessary logic for determining 
the pertubation points A and the orderly arrangement of the terms of 
the partials. The subprogram also performs the necessary 
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manipulations in calling the subroutine CRT which is used to deter- 
mine Cj^(jo)) in equation (A-16) . 

The input/output variables are defined in the following list. 
Note the compensator coefficients and related data enter the sub- 
program through a common block. 

Subprogram Variables 

Input Variables ; 

OMEGA C I) - A complex one dimensional array that contains the 

specified frequencies in ascending order for describing 
the system. 

NFREQ - An integer variable denoting the number of active 

margins to be improved. 

CT(I) - A one dimensional complex array represented by the 

compensated closed-loop frequency response C(s) in 
equation (A-19) . 

KPTS(I) - An integer array used as a pointer to denote the 
ftequency number of the margin investigated. 

TYPE (I) - A real array used to denote the type of margin being 

investigated. 

T(I,J,K) - A three dimensional complex array containing the 
transfer product of the compensation G(s) and the 
plant response P(s) for the specified frequencies, 
P(I,J,K) -A complex three dimensional array describing the open- 
loop frequency response of the plant. 

G(I,J,K) - A complex three dimensional array representing the 

compensation evaluated at the specified frequencies. 
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KSYM 

PFX(I,J), 

PFY(I,J) 

A,B, C^D , 
E,F,Z 


- An integer used for addressing the proper system 
arrays, 

- A two dimensional real array containing the partials 
of the numerator and denominator compensators respec- 
tively. 

- Parameter variables denoting dimension allocations. 


The following variables are used to describe the compensation 
polynomials; please refer to the Subprogram EVAL for a complete 
description; KIN,K0UT,ZA,ZB,ZC,ZD,ZE,PA,PB,PC,PD,PE,N1,N2,M1,M2, 


GAIN.KONT. 


The transient variables apply to the following auxiliary sub- 
program and are not affected directly in this routine: 

Subprogram CRT; C1,CI,W0RKI 


Output Variables : 

NPARC - An integer variable that represents the number of 

partials determined, 

PG(I,J) - A two dimensional real array representing the 3d/8w 


in (A-17) . 
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Subroutine PHASEM 

The subprogram PHASEM is used to detect and calculate the phase 
margins of an open-loop control system represented by a discrete 
frequency response. The open-loop frequency response is assumed to 
be given in terms of real and imaginary values. In particular given 

the ith frequency as f^, the corresponding real and imaginary parts 

of the frequency response are GR^ and The phase margins occur 

at the real zero crossings of the sequence: 

S. = 1.0 - jGR. + jGI.l^ . (A-24) 

X ' X x' 

Next the following sequence is formed: 


U. 

X 



(A-25) 


If TJ. < 0, then either S. or S. - is a zero or S. has made a 
zero crossing. Regardless of which condition has occurred, the 
frequency number of the phase margin is chosen as i or i-1 depending 
on the smaller magnitude of or The corresponding margin is 

calculated as 


STEM == |l.0 + GR^ + 

where k is either i or i-1 as mentioned above. 


(A-26) 


The following variables are defined for this program: 
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Subprogram Variables 

Input Variables ; 

OMEGA(I) - A complex one dimensional array that contains the 

specified frequencies in ascending order for describing 
the system. 

GTOTAL(I) - A complex one dimensional array containing the compen- 
sated open-loop frequency response for the Ith specified 
frequency point. 

KPOINT - The integer number of frequency points used to describe 
the open-loop frequency response of the system. 

FQMIN - The lowest frequency for which the particular margins 

are to be determined. 

KM - The integer used as a counter for the number of margins 

located.' 

Output Variables : 

KM - This is the number that designates the last margin 

found. 

STBM(I) - A one dimensional real array that contains the margin 
values corresponding to the frequency pointer KPTS (I) . 
These margins are measured in terms of distances from 
the (-1 + jO) point in the complex GH(jO)) plane. 
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Subroutine POLEV 

The purpose of this subprogram is to evaluated polynomials 
at specified frequency points. The frequency data is complex in 
nature. The routine is designed with an internal subfunction to 
avoid inaccuracies in raising a complex variable to a power. The 
input, output variables are defined as follows; 

Subprogram Variables 


Input Variables ; 

FW(I) - A one dimensional real array containing the poly- 
nomial coefficients in ascending order. 

K - An integer denoting the order of the polynomial. 

X - The complex variable of evaluation. 

Output . Variables : 

F - The complex resultant of the evaluation. 
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Subroutine SRMIMS 

The purpose of this subprogram is to calculate the maxima or 
minima of a discrete data open-loop frequency response with respect 
to a chosen point along the real axis. For example, assume GR^ and 
GI^ are the real and imaginary response terms corresponding to the 
ith frequency point. The following sequence is formed: 

= 1? + GR^ + jGI^i^ . (A-27) 

where P represents the negative point of investigation located on 
the real axis. Another sequence is generated as follows: 

V. = U. - U. - . (A-28) 

X X x-1 

.If * ^i 1 — ^i 1 ^ ^ ’ then the (i-1) frequency point 

corresponds to a relative maximum with respect to point P. Otherwise, 
if • ^i-1 less than or equal to zero and 0 » then the 

(i-1) frequency is a relative minimum with respect to P. 

The definitions of the variables are the same as those in the 
routine PHASEM with the following additional input variables: 

Subprogram Variables 

Input Variables : 

P - The negative value of the real axis point for which 

maxima or minima are to be located. 

N - An integer variable to determine whether the program 

is to determine maxima or minima. Maxima are investi- 
gated if N is unity; minima are found if N equals -1. 
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Subroutine XCHECK 

The XCHECK routine is not designed to check for damoine; ratio 
violations, but is necessary to assure that comOensator noles and 
zeros on the specified zeta boundaries do not result in a zeta 
violation upon incrementation. The technique employed for designing 
the compensation is to adjust the directional vector by zeroing the 
corresponding partial vector terms until the directions of movement 
of the above mentioned poles are within the defined sector. In this 
effort, the subsystem values of the zeta damping factor and the 
undamped natural frequency must be calculated and compared to the 
user specifications. 

Recall that the form of the compensation used by GIF is cascaded 
first and second order factors. In the case of first order factors 
if either of the coefficients is negative the program defaults with 
an error signal denoting the occurrence of a zeta violation, and, 
consequently, the run is temninated automatically; this can only 
occur on the first iteration because on succeeding iterations the 
subprogram, YCHECK, assures that a zeta violation cannot occur. The 
assumption here is that the user has erred in his initial selection 
of compensation. 

If either of the coefficients is zero, the result is a root at 
the origin or at infinity. In this case the avoidance of a zeta 
violation is assured by forcing the corresponding terms of the 
directional vector to be nonnegative. This is implemented by first 
checking the signs of the corresponding terms of the directional 
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vector; if these signs are negative, the corresponding terms of the 
partial vector are zeroed and the directional vector is recomputed. 
On each iteration the number of variable coefficients is reduced by 
the nxjmber of terms in the directional vector forced to zero in this 
manner. If the number of variable coefficients becomes less than 
the number of active margins the CIA will fail and the run will 
terminate automatically. 

In the case of second order factors the program calculates the 
and zeta as defined by the following second order factor: 


T(s) = + 2^0) s + ( 1 ) ^ . (A-29) 

n n 

This equation can. be related to the CIP form of second order factors 
as: 


T(s) = a^ + a^s + a^s^ 

Comparing (A-29) and (A-30) it is easily seen that 


(A-30) 


(A-31) 


C = (ap/(2o)^ • a^) . (A-32) 

The damping ratio ^ is checked to determine if a zeta violation 
has occurred. If such an occurrence is detected, as with first 
order factors the program is automatically terminated. 

If the value of C indicates roots of the second order factor 
on the Z boundaries, the next step is to determine whether these 
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roots will result in zeta violations upon incrementation of the 
compensator coefficients. This is accomplished by checking the sign 
of the change in 5 from incrementation. 

The change in zeta with respect to the change in the coeffic- 
ients of (A-30) is 



^ . 1 

.■75 Aa 4- — Aa, 

2^o^2 o a-2 ^ 



(A-33) 


where Aa^ , Aa^ , and Aag are the elements of the directional vector 
corresponding to aj, a^, and a^ . If A? is negative a zeta 
violation is inevitable. In order to avoid such an occurrence, 
associated terms of the partial vectors are zeroed and the directional 
vector is recommputed. This is continued until A^ is nonnegative. 
As with the first order factors the number of variable coefficients 
is reduced. 

The program variables are defined as follows: 


Subprogram Variables 

Input Variables t 

KIN - An Integer variable that denotes the number of 

inputs to the controller. 

KOUT - Integer variable denoting the number of controller 

outputs . 

ZA(I,J,K) - A real three dimensional array representing the 
constant terms of the first order factors. 

ZB(I,J,K) - A real three dimensional array containing the first 
order factor coefficients of s. 
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ZC(I,J,K) 

ZD(I,J,K) 

ZE(I,J,K) 

N1(I,J) 

N2(I,J) 

DV(I) 

ZETA 

KKK 

PG(I,J) 

KRE 

LPV 

NAM 

NRTRl 


- A real three dimensional array representing the 
constant terms of the second order factors. 

“ A three dimensional real array containing the s 
coefficients of the second order factors. 

- A real three dimensional array containing the s^ 
coefficients of the second order factors. 

- An integer array that denotes the number of cascaded 
first order factors for each subsystem. 

- An integer array that denotes the number of cascaded 
second order factors for each subsystem. 

- A real one dimensional array representing the 
directional vector ^ in equation (A-5) . 

- A real variable that denotes the desired minimum 
damping ratio. 

- An integer used to count the number of elements in 
the, directional vector. 

“ A two dimensional array containing the partial 
9^/9w in equation (A-17) . 

- A program control integer. 

“ An integer denoting the number of variables allowable 
for change. 

- An integer variable that represents the number of 
active, that is, unsatisfied margins. 

- An integer denoting whether first order factors are 

constrained to the left half GH(ja)) plane: if 1, the 

factors are unconstrained; otherwise, constrained. 
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NRTR2 - Same as NRTRl but applying to second order factors. 

A,B,C,D,E - Parameter variables used to dimension the arrays by 
the number of maximum allowable elements. 

Output Variables ; 

The following output variables are defined in the same manner 
as their respective input variables, but have been updated in the 
subprogram: PG(I,J), KKK, KRE, LPV. 
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Subroutine YCHECK 

The subprogram YCHECK is designed to detect compensation poles 
or zeros that have been forced outside the allowable camping ratio 
sector due to incrementation of the coefficients. If such is the 
case, the routine selects a maximum step size which will inhibit the 
zeta boundary violation while continuing to produce an improved 
solution. 

Recall that the CIP compensator data is described by transfer 
functions in cascaded first and second order factors. Similar to 
the XCHECK subprogram, YCHECK examines the first order compensator 
factors for possible right half plane roots. If such a root is 
found, the program reduces the step size until the root is marginally 
in the left half plane. 

Assuming that all the first order roots are now in the left 
half plane, the subprogram proceeds to- investigate the second order 
factors for possible zeta boundary violations. A t3T5ical second 
order factor of the form 

T(s) = ag + ajS + (A-34) 

yields the relation for the- undamped natural frequency 

% " , (A-35) 

and the zeta damping ratio as 


(aj)/(2w^ • a^) 


(A-36) 
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Upon comparison with the user specified damping ratio, if a zeta 
violation is incurred, the deincrementation of the second order 
coefficients continues iteratively until no violation occurs. From 
this increment a new step size is calculated and the routine returns 
to the main program to determine the compensator coefficients in 
accordance with this step size. 

The input and output variables are defined in the following 

list. 


Subprogram Variables 

Input Variables : 

KIN - An integer variable that denotes the number of inputs 

to the controller. 

KOUT - Integer variable denoting the number of controller 

outputs . 

ZA(I,J,K) -> A real three dimensional array representing the con- 
stant terms of the first order factors. 

ZB(I,J,K) - A real three dimensional array containing the first 
order factor coefficients of s. 

ZC(I,J,K) - A real three dimensional array representing the con- 
stant terms of the second order factors; i.e., a^ 
in (A-34) . 

ZD(I,J,K) - A real three dimensional array containing the s coef- 


ficients of the second order factors; i.e., in 
(A-34) . 
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ZE(I,J,K) 

N1(I,J) 

N2(I,J) 

DV(I) 

ZETA 

KKK 

STEP 

PMG 

KRE 

A,B,C,D,E 

NRTRl, 

1S1RTR2 


A real three dimensional array containing the 
coefficients of the second order factors; i.e., 
in (A-34) . 

An integer array that denotes the number of cascaded 
first order factors for each subsystem. 

An integer . array denoting the number of cascaded 
second order factors for each subsystem. 

A real one dimensional array representing the 
directional vector jd in equation (A-5). 

A real variable that denotes the desired minimtim 
damping ratio. 

An integer used to count the number of elements in 
the directional vector. 

A real variable denoting the step size . 

A real variable denoting the magnitude of the partial 
vector. 

A program control integer. 

Parameter variables used to dimension the arrays by 
the maximum number of elements allowable. 

An integer denoting whether first or second order 
factors, respectively, are constrained to the left half 
plane; if 1, the factors are unconstrained; otherwise, 
constrained. 


Output Variables 

The following variables are defined in the same manner as their 
respective input variables, but have been updated in the subprogram: 


STEP, KKK, KRE, LPV. 



APPENDIX B 


FORTRAN IV LISTING OF THE COMPENSATOR 
IMPROVEMENT PROGRAM 

This appendix contains a complete Fortran version of the 
Compensator Improvement Program for Multi-variable Control Systems. 
The program is completely self-contained, i.e., it requires no 
system library, etc. The necessary data input is explained in the 
comment statements preceding the main program. The subprograms are 
listed in alphabetical sequence. For information regarding sub- 
program theory or variables, refer to the respective synopsis of 
Appendix A. 



Ill 
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11 + NG + NP + NS ASF(L),ASF(L),L=UNA 8G10.5 


• 

12+NG+NP+NS+NA FlrF2/ . . ,/FlO 8G10.3 

OMEGACK) ,G(I,J,K) ,J=1,KIN, 

, I=1,K0UT aG20.5 

GAINCl, J) ,N1 (I,J] ,N2C1,J) , 

. Ml (1/ J) »M2(I f J) ,KONT(I, J) GIO.5,515 

. ZA(I,J,L)»ZB(I,J,U,L = i/Nl 7G10.5 

ZCCI, J,L) /ZDCI/J/U ,ZE(I,J,LD , 

. L=1,N2 7G10.5 

. PACI, J,L)/PB(1,J,L),L=1,M1 7G10.5 

PCU,J,L),POa,J,L)»PEa, J,L)f 
. L=1,M2 7G10.5 

DEFINITIONS OF I/O VARIABLES 

ITERATION CONTROL DATA 

MODE -DETERMINES WHETHER THE PROGRAM IS TO OPERATE IN THE 

SUM improvement frouemcy mode (sifr) or the total 
improvement frequency mode CTIFRI. if data CARO 

BLANK THE PROGRAM AUTOMATICALLY DEFAULTS TO THE 
TIFR MODE. IF SIFR IS DESIRED THEN MODE"SIFR, 
NZEROl-DETERMIMES VnHETHER 1-ST ORDER ZERO FACTORS ARE 

constrained to L.H.P. (if .EQ. TO 1, UNCONSTRAINED 
OTHERiMSE, CONSTRAINED). 

NZER02-SAME. AS NZEROl, EXCEPT FOR 2-N‘l) ORDER FACTORS. 
NPOLEl-SAME AS NZEROl, EXCEPT FOR 1-ST ORDER POLE FACTORS. 
NPOLE2-SAME AS NZEROl EXCEPT FOR 2-ND ORDER POLE FACTORS. 
ZETAZ- MINIMUM DAMPING RATIOS FOR COMPENSATOR ZEROS; APPLIES 
ONLY IF NZERO .NE.l. 

KEY = 0 , NYQUIST SUBPROGRAM NOT CALLED 
ZETAP - SAME AS ZETAZ APPLIED TO POLES 
KSTART -STAH.TING ITERATION NO. 

KWUIT - STOPPING ITERATION NUMBER 

KPOINT -NO. OF POINTS FROM OPEN LOOP FREQ. RESPONSE USED 

IF KPOINT=0, READ A FREQUENCY, THEN EACH CHANNEL'S 
STPMAX -MAXIMUM CHANGE TO BE MADE IN COMPENSATOR COEFFICIENTS 
SMALLEST COMPENSATOR COEFFICIENT OF THE INITIAL 
COMPENSATOR) 

STPMIN - MINIMUM STEP SIZE DESIGNATOR 

PINACT -LARGEST DIFFERENCE BETWEEN A CONSTRAINT AND ITS 

DESIRED VALUE IN GOING FROM INACTIVITY TO ACTIVITY 
KPRINT - MO. OF ITERATIONS SKIPPED BETWEEN PRINTING OF IMFOR, 

KEYOUT=l , REQUESTS PRINTOUT 
KEYOUTsO , NO PRINTOUT 
KEYOUTCl) complete ITERATION OUTPUT 
KEY0UT(2) - OUTPUT COMPENSATOR INFORMATION 
KEY0UTC3) - OUTPUT FREQUENCY RESPONSE 
KEYOUT(A) - OUTPUT MARGIN SPECIFICATIONS 
NONCTRCI) - NUMBER OF POLES ON CONTOUR 
NOLRHPCI) - NUMBER OF POLES INSIDE CONTOUR 
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VARIABLES FOR MARGIN RADII SPECIFICATIONS 

NUMBER OF MARGIN RADII TO BE SPECIFIED; 

NG - NO. OF GAIN MARGIN RADII SPECIFIED 
NP - NO, OF PHASE MARGIN RADII SPECIFIED 
NS - NO. OF STABILITY MARGIN RADII SPECIFIED 
NA - NO. OF ATTENUATION MARGIN RADII SPECIFIED 


VARIABLES FOR GAIN MARGIN RAD.II DESIGNATIONS 
IF FREQ .GT. GMFCI) BUT .LT. GMFCI-rl) DESIRED MARGIN = GMR(I) 

IF FREQ .GT. GMF(I+n BUT .LT. GMF(I+2) DESIRED MARGIN s GMRCI+I3 
ETC, 

IF FREQ .GT, GMF(NF) DESIRED MARGIN = GMR(NF) 

VARIABLES FOR PHASE MARGIN RADII DESIGNATIONS 
IF FREQ ,GT. PMFCI) BUT .LT. PMF(I+n DESIRED MARGIN s PMR(I] | 

IF FREQ ,GT.PMFCI + n BUT ,LT, PMFCI+2) DESIRED MARGJN s PMR(I + n 

ETC. 

IF FREQ .GT, PMF(NF) DESIRED MARGIN = PMR(NF) i 

VARIABLES FOR STABILITY MARGIN RADII DESIGNATIONS 
IF FREQ .GT, SMF(I) BUT .LT. SMF(I+1) DESIRED MARGIN = SMR(I) ^ 

IF FREQ .GT. SMFCI + 11 BUT .LT. SMFCI+2) DESIRED MARGIN = SMRU + l] 

ETC, 

IF FREQ .GT, SMF(NF) DESIRED MARGIN = SMR(NF) 

VARIABLES FOR ATTENUATION MARGIN RADII DESIGNATIONS 
IF FREQ .FT. ASFd) BUT .LT. ASFf.l + l) DESIRED MARGIN = ASRCIJ 
IFCFREQ .GT. ASF(l+n BUT .LT. ASFCI+2) DESIRED MARGIN = ASRCI) 
ETC-. 

IF FREQ .GT. ASF{NA) DESIRED MARGIN = ASRCNA) 


FI ; F2 - FREQUENCIES, BETWEEN WHICH G.M.'S ARE FOUND 
F3 ; FA - FREQUENCIES BETWEEN WHICH P.M.'S ARE FOUND 
F5 ! F6 - FREQUENCIES BETWEEN WHICH S.M.'S ARE FOUND 
F7 ! F8 - FREQUENCIES BETWEEN WHICH A.M.'S ARE FOUND 


VARIABLES FOR PLANT DYNAMICS 

KIN IS NO, OF CONTROL INPUTS' 

KOUT IS NO. OF OUTPUTS 

OMEGACI) - ITH FREQ. C ASSUMED TO BE I-N HZ.) 

G(I,J/K) - PLANT DYNAMICS WITH INDEX F CK IN , K OUT , KPOl NT ) 

DESCRIPTION OF COMPENSATION 

GAINCn-DENOTES INITIAL D. C, GAIN VALUE FOR I"TH CHANNEL 
KONT(I)-D.C, DESIGNATOR FOR I-TH CHANNEL 

KONT(I)=l GAIN ALLOWED TO VARY 

K0NT(I)=2 GAIN NOT ALLOWED TO VARY 




N1(I,J3 - MO. OF iST order MU^iEWATOR COEFFS. OF COMPENSATOR C1»J) 

- NO. OF 2ND ORDER NUMERATOR COEFFS, OF COMPENSATOR (I/J) 

- NO. OF 1ST ORDER DENOMINA. COEFFS. OF COMPENSATOR {IrJJ 

M2(I,J) - NO. OF 2ND ORDER DENQMINA. COEFFS, OF COMPENSATOR (I»J) 

KTH CASCADE COMPENSATOR ( I , J , K ) COEFFICEMTS .^JHERE I = KIN, J = KOUT 
FIRST ORDER NUMERATOR FACTOR COEFFICIENTS: (ZA + ZB*S). 

ZA(I,J,K) - 1ST ORDER NUN. FACTOR CONSTANT OF KTH CASCADE COMP, 
ZB(I,J,K) - 1ST ORDER NUMERATOR FACTOR COEFF IC lENT C ZB*S) OF KTH 
SECOND ORDER NUMERATOR FACTOR COEFFICIENTS: (ZC + ZO*S + ZE*S**2) 

ZC(I,J,K) - 2ND ORDER MUM. FACTOR CONSTANT OF KTH CASCADE COMP. 
ZO(I,J,K) - COEFFICIENT ZD*S OF KTH CASCADE COMPENSATOR 
ZECIrJ,Kl - coefficient ZE*S**2 Of- KTH CASCADE COMPENSATOR 
FIRST ORDER DENOMINATOR FACTOR COEFFICIENTS.: CPA + PB*S) 

PACI,J,K) “ 1ST ORDER DENOMINATOR FACTOR CONSTANT 
PB(I,J,K) - 1ST ORDER DENOMINATOR FACTOR COEFF, (PB*S) OF KTH ' 
SECOND ORDER DENOMINATOR FACTOR COEFFI C lENTS ! C PC+PD*S+PE*S**2 ] 

PCCIrJ>K5 - 2ND ORDER DENOM. FACTOR CONSTANT OF KTH CASC, COMP. 
PD(I,J,K} - coefficient PD*S OF KTH CASCADE COMPENSATOR 
P£CI,J,K) - COEFFICIENT PE*S**2 OF KTH CASCADE COMPENSATOR 
INTEGER TYPEf TYACT, ACTIVE, ACTDES,XXP,XXG/XXS,XXA 
COMPLEX G(8,A,'D5 #GC(A,B,D),T(A',A,0) ,OM£GACD).,CT(0,A) ,.PCG,C1 (A,A),C 
1I(A,A),W0RKI(A,A2) 

OIM£NSION K0NTCA.,B), GAINCA/B), N1(A,B), M2(A,8), M1(A,B1, M2(A,BJ 

1, ZA(A,B,C), ZeCA,B,C), ZC(A,8,C), ZD(A,B,C), ZE(A,8,C), PACA,B,C) 

2, PBCA,B,C), PC(A,B,C), PD(A,B,C), P£(A,B,CD, GMFCH), GMR(H), PMF ( 
3H), PMR(H), SMF(H), SMRCHJ, ASFCH), ASR(H), NMCA,, KOLDCOl, PG(E,F 
«), DVCF), ^EIGHTCE), KMINCAl, SKL(E,A1, JYACT(E), RO(E,AT, TYPECE, 
5A), ACT1V£CE,A1, ST8M(E,A), KPT.S(E,A], NIT(A), KINACT(E,A), NMLCA) 
6, KACTCE,A), IDC15), ACTDE5(2), PFX(E,E3, PFYCE/E3, WORKR(E,E2), K 
7NErt(0), NONCTRCAI, NOLRHP(A), KEY0UT(4) 

COMMON 1TER,KSTART 

DATA IBLANK,XXP,XXG,XXS,XXA /4H , 1 HP , 1 HG , 1 HS , 1 HA/ ACT0E3Cn,ACT 

1D£S(23 , ITIFR /3HYES, 2HNO, UHTIFR/ 

KPTMAX=0 

DATA INPUT BLOCK -T^ 

ITERATION CONTROL DATA 

READ (5,10) MODE,NZER01,NZER02,NPOLE1,NPOLE2,KEY,Z£TAZ,ZETAP,ID 
10 FORMAT (Aa,6X,5I5,2G10,5/15A4) 

VYRITE (6,20) ID 
20 format (IHl , lOX, 15A4) 

WRITE (6,30) N’ZEROl ,NZEP02,NPOLE1,NPOLE2 
30 FORMAT C'O^SX, 'NZER01 k',I2,' NZF.R02= ' , 12, ' NP0LE1 = ", 12, ' NP0LE2 
1=',12) 

IF (MODE.EQ. IBLA'NK) MODE=ITIFR 
WRITE (6,40) mode 

40 FORMAT {'0',5X,'THE PROGRAM IS IN THE ',A4,' MODE') 

READ (5-,50) KSTART,KaUIT,KPOINT,KPRINT,KIM,KOUT,STPMAX,STPMIN,PINA 
ICT 

so format (6I5/3F10.5) 

READ (5,60) (KEYOUT(T),I=l,4) 

READ (5,, 60) (NONCTR(l) ,N0LRHP(I)., 1 = 1 ,KIN) 

FORMAT (615) 


60 




115 



DESIGN SPECIFICATION DATA , 

READ (5r60) NG/MP,NS,NA,NV 
WRITE C6,80D KST ART , KQUI T » KPOI NT , KPRINT 
READ MARGIN FREQUENCIES IN ASCENDING ORDER 
READ (5,70) (GMF(K),GMK(K),K=1,NG> 

READ (5,70) (PMF(l) ,PMR(L) ,L=1 ,NP) 
read (5,70) (SMF(1< ) ,SMR(K) ,K = 1 ,NS) 

READ (5,70) (ASF(L),ASR(L),L=1,NA) 

READ (5,70) FI ,F2,F5,F4,F5,F6,F7,F8,F9,F10 
FORMAT (BG10,3) 

FORMAT ('O', 5X, 'START ITER = ',IA,' STOP ITER = ',IU,2X,'N0. OF 
IFREQ. POINTS = ' , IS , 2X , 'PR IM INCREMENT = ',15) 

WRITE (6,90) KINjKOUT 

FORMAT ('O', 5X, 'NUMBER OF CONTROL INPUT CHANNELS= ' , 15 , 5X , ' NUMBER 
lOUTPUT CHANNELS = ',15,/) 

WRITE (6,100) stpmax,stpmin,zetaz,zetap,pinact 

FORMAT ( '0 ',5X, 'MAXIMUM DESIGNATED STEP SIZE =', F 1 0 . 5 '6X , 'MINIMUM 
IDESIGNATED STEP SIZE = ', G 1 0 .5/6X , 'MINIMUM COMPENSATOR ZEROS ZETA 
2= ' ,G10.5/6X, 'MINIMUM COMPENSATOR POLES ZETA s ' , G1 0 , 5/ 6X , ' AMT . A 
30VE SPECS, TO BE KEPT IN ACT , = ' , F 1 0 , 5 ) 

DO no 1 = 1, KIN 

WRITE (6,120) I,NONCTR(I),NOLRHP(I) 

format ('0',5X,'F0R system NO.', 13,'; NONCTR=' , 13, ' AND NOLRHP=', 
13) 

WRITE (6,130) NG,NP,NS,NA 

FORMAT ( '0',5X, 'NG='I2, ' ,NP=',I2,' ,NS=',I2,' ,NA=',I2) 

WRITE (6,ia0) 

FORMAT ( '0'5X, 'DESIRED MARGIN RADII DESIGNA-T IONS ', ///6X , 'DESIRED 
I'GAIN MARGIN DESIGN SPECIFICATIONS') 

WRITE (6,150) (GHF(K),GMR(K),K=1,NG) 

FORMAT (/,15X,'IF FREQUENCY . GE . ' , F 1 0 . 5 , 5X , ' DESI RED MARGIN IS', FI 
1.5) 

WRITE (6,160) F1,F2 

FORMAT C' ',5X,'GAIN MARGINS ARE DETERMINED BETWEEN THE FREOUENCI 
IS OF',F10.5,2X, 'ANO',F10.5) 

WRITE (6,170) 

FORMAT ('0',5X, 'DFSIRED PHASE MARGIN DESIGN SPECIFICATIONS') 

WRITE (6,150) (PMF(L),PMR(L),L=1,NP) 

WRITE (6,180) F3,F« 

FORMAT (' ',5X, 'PHASE MARGINS ARE DETERMINED BETWEEN THE FREOUENCI 
lES OF',F10.5,2X, 'ANO',F10.5) 

WRITE (6,190) 

format ( 'O', 5X, 'DESIRED STABILITY MARGIN DESIGN SPECIFICATIONS') 
WRITE (6,150) (SMF(K) ,SMR(K) ,K=1 ,NS) 

WRITE (6,200) F5,F6 

format (' ',5X, 'STABILITY MARGINS ARE FOUND BETWEEN TH£ FREQUENCIE 
IS OF',F10.5,2X, 'AND',F10.5) 

WRITE (6,210) 

FORMAT ('O', 5X, 'DESIRED ATTENUATION MARGIN DESIGN SPECIFICATIONS') 
WRITE (6,150) (ASF(K),A5R(K),K=1,NA) 
vJPITE (6,220) F7,F8 

FORMAT (' ',5X, 'ATTENUATION MARGINS ARE FOUND BETWEEN THE FREQUENC 
IIES OF',F10.5,2X,'ANO',F10.5) 

READ DATA ON O.L. SYSTEM 

I,J,K ALWAYS DENOTE COUNTER K IN, KOUT , KOATA POINT RESPECTIVELY 


,NS=',I2,' ,NA=',I2) 
DESIGNATIONS', ///6X, 'DESIRED 


SPECIFICATIONS') 
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WRITE (6,230) 

30 FORMAT (5X,'0PEN LOOP SYSTEM INPUT FREQUENCY RESPONSE: ',/,31X,'CO 
IMPLEX OMEGA', 13X, 'complex G(JW) '/,5X, 'DATA ', 3X , 'CHANNEL ( I , J) ' ,5X , ' 
2REALCW) ',5X, 'IMAG(W) ',5X, 'REAL(G) ' ,5X, ' IMAG{G) ',/) 

00 260 K=1,KP0INT 

READ (5,200) OMEGA(K), C(G(I,J,K),J=I,KIN),I=1,K0UT) 

00 format (OG20.5) 

00 260 I=1,K0UT 
DO 260 J=1,KIN 

aRITE (6,250) K,I,J,OMEGA(K),G(I,J,K) 

250 FORMAT (5X,I3,2(5X,I2),2(5X,2G10,«)) 

260 CONTINUE 

C COMPENSATOR INPUT DATA . . , 

DO 320 1=1, KIN 
DO 320 J=l,KOUT 

READ (5,270) GAINCI,J),NUI, J) ,N2(I,J) ,M1 ( I , J) ,M2 (I , J) , KONT C I , J ) 
270 FORMAT (IGIO.5,515) 

C . read COEFFICIENTS IN ASCENDING POWERS OF S, IN 1ST AND 2ND ORDER 
NC=Ml(I,J) 

IF (NC.EQ.O) GO TO 280 

READ (5,310) (ZA(I,J,L),ZB(I,J,L),L=1,NC) 

280 NC=N2(I,J) 

IF (NC.EQ.O) GO TO 290 

READ (5,310) (2Ca,J,L),ZD(I,J,L),ZE(I,J,L),L=l,NC) 

290 MC=M1(1,J) 

IF (MC.EO.O) GO TO 300 

READ (5,310) (PA(1,J,L) ,PB(I,J,L) ,L=1,MC) 

300 MC=M2(I,J) 

IF (MC.EQ.O) GO TO 320 

READ (5,310) (PCa,J,L),PO(I,J,L),PE(I,J,L),L = l,MC) 

310 FORMAT (7G10.5) 

320 CONTINUE 

C DATA INPUT COMPLETED 

|C determine no. of INDEPENDENT PARAMETERS 

KN0T=0 
LNOT=0 
KVARYsO 

DO 330 1=1, KIN 
DO 330 J=1 ,KOUT 
KNOT=KNOT+2*N1 (I, J)+3*N2(I, J) 

LN0T=LN0T+2*M1 (I, J)+3*M2(I, J) 

KVARY=KVARY+N1 (I,J)+N2(I,J)*2 
KVARY=KVARY+M1 (I, J)+M2( I , J)*2 
330 IF (K0NT(I,J) .EO.l) K VARY=KVAR Y+ 1 

nparc=knot+lnot 

iter=kstart 

step=stpmax 

STPOLO=STPMAX 

kpold=kpoint 

DO 3il0 I = 1,KP0LD 
3«0 K0LDCI)=I 
KPR=-1 

DATA RAD2 /lU. 591559/ 

DO 360 1=1, NP 

IF ((PMR(I).GE.O.).AND.(PMR(I),LE.180.)) go to 360 
WRITE (6,350) 

350 FORMAT ( '0 ', 5X ,.'**** HEY DUMMY YOU HAVE MADE A MISTAKE" ON THE 

IPHASE margin specifications ****') 
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STOP 

360 PMR (na2,*SIN(PMR(I)/RAD2D 

FNY0T = AMAX1 (GMF CNG),PMF CNP) ,SMFCNS) ) 

370 KRESET=0 

IF { (ITER. EQ.KSTART3 .OR. (KPOINT.lt, 1 .5*KSTAND) ) GO TO 3SO 
CALL DELETE ( KPOINT ,K IN, KOUT, OMEGA, G, NIT , KINACT , NML, K ACT, ITER, KPOl 
10,K0L0,KNEW,A,B,C,D,E) 

KRESET=1 

C CALCULATION OF COMPENSATED FREQUENCY RESPONSE ... 

380 DO 390 K=l,KPOINT 

390 CALL EVAL (OMEGA (K ) , GC , G, T , K, K IN , KOUT , Z A , ZB , ZC , ZD, Z£ , PA , P8 , PC , PO , I 
IE,N1 ,N2,M1 ,M2,DAIN,K0NT,A,8,C,D) 

IF (KPR.EQ.KPRIMT) KpR=0 
KPR=KPR+1 

IF (ITER.NE.KSTART) GO TO 4J00 

IF (KEYOUT(2) ,EQ. 1 ) CALL OUTPT ( 0 , KS YM, KPOINT , CT ( 1 , KS YMD , OMEGA , I Tl 
1R,N,KPTS(1 ,KSYM) ,KMIN,STBM(1,KSYM) ,RQ,TYPE(1 ,KSYM) , ACTIVE,KIN, 

2 K0UT,ZA,Z8,ZC,ZD,ZE,PA,PB,PC,PD,PE,N1,N2,M1 ,M2,GAIN.K0NT,A,B,C,! 
30LDI 

IF (ITER.GT.KQUIT) go to 1020 

C determination - OPEN LOOP FREQUENCY RESPONSE CT OF K3YM .... 
400 NAM=0 
NSUMsO 
MA0 = 0 
LPVsKVARY 
ADD=1 .E-15 
DO 750 KSYM=1,KIN 
KPLASTsKPOINT 
DO 410 K=1,KP0INT 

CALL CRT (l,K,KSYM,T(l,i ,K3 ,KIN,KOUT,CT(l ,KSYM) ,G(l, I ,K) ,PCG,I1 , J 
1 ,C1 ,CI,W0RKI, A,B,D) 

IF (KEYOUT{33.EQ.n CALL OUTPT ( 1 , KSYM, KPO INT , CT ( I , KSYM) , OMEGA , IT 
1R,N,KPTS(1 ,KSYM) ,KMIN,STBMCi,KSYM) ,RQ,TYPE ( 1 , KSYM) , AC Tl VE , K IN , 

2 K0UT,ZA,Z8,ZC,Z0,ZE,PA,PB,PC,PD,PE,N1 ,N2,M1 ,M2,GA1N,K0NT,A,B,C, 
30LD) 

10 CONTINUE 

DETERMINATION OF MARGINS . , . 

20 NM(KSYM)=0 

DETERMINATION OF GAIN MARGINS BETWEEN FI AND F2; 

CALL GAINMG (CT ( 1 , KS YM3 , KPOINT , NM (KSYM) , F 1 , F2 , KPTS ( 1 , KSYM3 , STBM ( 1 
IKSYM) .OMEGA) 

KPM=NM(KSYM)+1 

NGMS=NM(KSYM) 

IF (NMCKSYM) ,EQ.0) go TO 440 

call adopts (KP01NT,KIN,KOUT,l,NM(i) , STBM( 1 , K3 YM) , KPTS ( 1 , I ) , CT ( 1 , 1 
I ) ,G,GC,T,OMEGA,KGOHAK,KPTMAY,NIT,KINACT,N-'^L,KA' T(1 , 1 ) ,KPOLD,KOLD,K 
2SYM,C1 ,CI ,WORKI,ZA,ZB,ZC,ZO,2E,PA,PB,PC,PO,P£,N1 ,N2,M1 ,M2,GAIN,K0N 
3T,A,B,C,0,E) 

IF (KPOIMT.GT.KPTMAX) go to 450 
IF (KG08AK.EQ. 1 ) GO TO 420 
N=NM(KSYM) 

C SETTING DESIRED STABILITY RADII OF GM'S; 

DO 430 1=1, N 

TYPECI,KSYM)=XXG 

KWHICH=KPTS(I,KSYM) 

FR£HZ=AIMAG{0MEGACKWHICH) ) 
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DO «30 L=irNG 

IF CFREHZ.GE.GMFCLD J RQ ( I , KSYM) =GMR (L ) 

CONTINUE 

NM(KSYN)=NGMS 

•determination of phase margins between F3 A^fD F4; 

CALL PHASEM CCT ( 1 , KS YM ) , KPOIM , NM (KS YM ) , F3 » F4f KPTS ( 1 » KSYM) , STBM U , 
IKSYM) , OMEGA) 

N=NM(KSYM) 

IF (n.lt.kpm) go to 4R0 

CALL ADOPTS (KP0INT,KIN,K0UT,KPM,NMCn,ST8Mn ,KSYM) rKPTSd / 1),CT(1 
1 , 1 ) ,G,GC,T,0MEGA,KG08AK,kPTMAX,NIT,KINACT,NML/ KACTCd 1 ) ,KPOLO,KOLD 
2,KSY.M,C1 / CI/rtORKI,ZA,ZB,ZC,ZO,ZE,PA,P8,PC,PD,PE,Nl',N2,Ml,Ma,GAIN,K 

30NT,A,8,C,D,E) 

IF CKPOINT.LE.KPTMAX) GO TO 470 
WRITE (6/460) KPTMAX 

format C '0 ',5X,3C1H*) /2X/ 'TERMINATION REASON: NO, OF FREQ. "POIN 

ITS HAS EXCEEDE0'/I5/2X,3C1H*)) 

GO TO 900 

IF CKGOBAK.EQ, 1 ) GO TO 440 
N=NM(KSYM) 

SETTING DESIRED STABILITY RADII OF P.M.'S; 

DO 480 I=KPM,N 
TYP£CI,KSYM)=XXP 
KWHICH=KPTS(If KSYM) 

FREH2=AIMAC(0MEGA(KWHICH) ) 

DO 480 L=1/NP 

IF (FREHZ.GE.PMF(D) RQ ( I ,KSYM) =PMR CL) 

CONTINUE 

KPM=N+1 

CONTINUE 

KLASTsNM(KSYM) 

NM(KSYM)=KLAST 

KST8M=KPM 

determination of STABILITY MARGINS 

call SRMINS (C-T(1,KSYM) /KPOINT/NMCKSYM) , I . / 1 , F5 / F6 , KPTS (1 / KSYM) , ST 
iaM(l,KSYH)/OMEGA) 

N=NM(KSYM) 

IF (N.LT.KPM) GO TO 590 

CALL ADOPTS CKPOINT , KIN , KOUT, KPM, NM (3 ), STBM Cl , KSYM) , KPTS C 1 / 1 ) ,CT(1 
1 / 1 ) ,G,GC/T,0!1EGA,KG08AK,KPTMAX,N1T,KINACT,NML/KACT(1 / 1 ) ,KPOLD,KOLD 
2,KSYM,C1 ,CI,WORKI,2A,Z8,ZC,ZO,ZE,PA,P8,PC,PDrPE,Nl ,N2,M1 ,M2,GAIN/K 
iONT,A,B,C/D,E) 

IF (KPOINT.GT.KPTMAX) go to 450 
IF (KGOBAK.ECJ. 1) GO TO 500 
IF (I TER.EQ.KSTART) go to 510 
IF (KEY.EQ.O) 60 TO 530 

CALL NYQIST (KPOINT/CT (1 , KSYH) , NOLRHP ( KSYM) , NONCTR (K S YM) / NC I RL / NZ ; 
IFMYOr, OMEGA) 

IF (NZ.EQ.O) GO TO 530 
IF (ITER.NE.KSTART) go to 760 
aRITE (6,520) KSYM 

FORMAT ( '0' ,5X, 'KITH THE INITIAL COMPENSATION FOR SYSTEM NO. ',13,1 
IX/'HAS CLOSED LOOP POLES INSIDE THE CONTOUR IN THE" PHASE STACilLI 
2ZATI0N REGION') 

CONTINUE 

IF (KPR-H .NE.KPRINT) GO TO 530 
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C SETTING DESIRED STABILITY MARGINS 

N=NM(KSYM) 

00 5«0 I=KPM,N 
TYPECI,KSYM)sXXS 
KvmiCH=KPTSCX ,KSYM3 
FR£HZ = AIMAG00MEGACK1\HICH3 ) 

DO 540 L=1,NS 

IF (FREHZ.GE.SMF(U ) RQ (1 ,KSYM)=SMR CLD 
540 CONTINUE 

C CHECK TO SEE IF ANY P.M, 'S, G.M. 'S# OR S.M^'S EQUAL 
C IF THERE RESULTS SOME THAT ARE EQUAL ONLY THE FIRST IS RETAINED. 

IF (KSTBM.LE.n GO TO 580 
DO 570 LB=l,N 
I NSG=LB+1 

!550 CONTINUE 

IF (NSG.GT.MM(KSYM) 3 GO TO 580 
DO 570 I = NSG/.N 

IF (KPTSCL8',KSYM3 .NE.KPTS(I,KSYM33 GO TO 570 
N=N-l 

DO 560 L=L8,N 

KPTS(L.KSYM3=KPTS(L+I/KSYM3 
ST8M(L,KSYM3=STBM(L+1 ,KSYM3 
RQ(L,KSYM)=RQ(Ltl»KSYM3 
560 TYPE(L,KSYM)=TYPE(L+1,KSYM3 
GO TO 550 
570 CONTINUE 
580 CONTINUE 
KPM=N+1 

590 CONTINUE 

KPAODsKPOINT-KPLAST 

IF (KPADD.NE,03 WRITE (6,6003 KPADO , I TER , KPO INT 
600 FORMAT C 1 HO , 5X , 1 3 , 1 X , 54HP0INTS WERE ADDED ON ITERATION NO.,I4,/10 
1,26HN0. OF FREQ. POINTS IS NOW, 143 
KMIN(KSYM3=N 

IF ( (ITER.EO.KSTART3 .OR. CKRESET.EQ.l 3 3 KSTANO=KpOINT 
C DETERMINATION OF ATTENUATION MARGINS 

CALL SRMINS (CT(1 ,KSYM),KP0INT,NM(KSYM3 ,0. ,-l,F7,F8,KPTS(l ,KSYM3 , 
ITBM(1,KSYM3,0MEGA3 

C SETTING DESIRED STABILITY MARGINS . , . , 

N=NMCKSYM3 

IF CN.LT.KPM3 GO .TO 620 
00 610 ,I=KPy,N 
TYP£(I,KSYM3=XXA 
KWHICH=KPTS(I,KSYM3 
FREHZ=AIMAG (OMEGACKWHICH33 
00 610 L=1,NA 

IF (FREHZ.GE.ASFCL33 RQ ( I , KSYM3 =ASR (L 3 
610 CONTINUE 
620 CONTINUE 

C CHECKING MODE requirements . . ..... 

IF (ITEP.EQ.KSTART3 GO TO 750 
P0RM=1 , 

N=NM(KSYM3 
N1=NIT(KSYM3 
NL=NML(KSYM3 
IF (N.E0.03 GO TO 750 
KCHMIN=MIN0 (N,NL+NI3 
KKK = 0 




00 660 Isl/N 
DO 630 J=1,NL 
DO 630 K=-2,2 

630 IF CKPTS(I,KSYM) ,EQ.KACTCJ,KSYM)+K) GO TO 650 
DO 640 J=1,NI 
DO 640 K=-2,2 

640 IF (KPTSCI,KSYM),EQ.KINACt(J/KSYM)+K) GO TO 650 
GO TO 660 
650 KKKsKKK+1 
660 CONTINUE 

IF (KKK.LT.KCHMIN) go TO 760 

nsum=nsum+n 

DO 740 1 = 1 , N 

IF (I.GT.KMINCKSYMI ) P0RM=-1, 

IF (PORM*(ST8M(I,KSYM)-RQn,KSYM)).GE.O.) GO TO 740 
IF (NIT(KSYM) .EQ.OO GO TO 700 
DO 690 J=1,NI 
DO 670 K=-2,2 

670 IF (KPTSCI ,KSYM) .EQ.KINACT (J,KSYM)+K) GO TO 680 
GO TO 690 

680 IF CPORM*CRO(I,KSYM)-STBMn,KSYMn .GT.PINACT) GO TO 760 
690 CONTINUE 

700 IF (NL.EQ.O) GO TO 740 
DO 710 J=1,NL 
DO 710 K=-2,2 

710 IF (KPTSCI, KSYM),EQ.KACT(J,KSYM)+K) GO TO 720 
MAD=MAO+l 
GO TO 740 
720 CONTINUE 

IF (MODE.NE.ITIFHD go to 730 

IF (PORK*(STBMa,KSYM)-SML{J,KSYM3).LT.-l,E-06) GO TO 760 
730 CONTINUE 

AD0=ADD+P0RM*(STBM(I,K3YM)*SML(J,KSYM) ) 

740 CONTINUE 
750 CONTINUE 

IF (ITER.EQ.KSTART) GO TO 800 
IF (MAD.EO.NSUM) ADD=.l. 

IF CADD.G7.0.) GO TO 770 
760 STEP=-ARSCSTEP}/2, 

IF CAas(STEP) .LT.STPMIN) GO TO 780 

ITER=ITER-l 

KPR=KPR-1 

GO TO 980 

770 STEP = l,414i6-*ABSCSTP0L'D) 

IF (ARS(STEP) .GT.STPMAX) STEP = STPMAX 
GO TO 800 

C OUTPUT CONTROL 

730 WRITE (6,790) STPMIN; 

790 FORMAT C '0 ',5X , '*,*** TERMINATION - STEP SIZE LESS THAN ',G15.5,2X 

1 '’****' ) 

GO TO 900 

800 NAM1=0 

IF ( CKPR.EO.KPRINT) .OR. ( I TER . EQ .KST ART) , OR , ( ITER.EQ.KQUIT) ) WRITE 
1(6,810) STEP 

810 format ('O', 25X, 'PRESENT STEP SIZE =',G15.5) 

DO 840 KSYM=i,KIN 
N=NM(KSYM) 

KSM=KMIN(HSYM) 
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OETEf^MINING ACTIVE MARGINS 

NIT CKSYM)=0 
K = 0 

DO 830 I=lfN 
P0RM=1 , 

IF (I.GT.KMINCKSYM)) porm=-i,. 

IF (PORM*STBM(I/KSYM) ,LT.PORM^tRQ(l,KSYM)) NAM1 = 1 
IF (P0RM*STBM(I,K3YM) ,GT.P0RM*HQCI#KSYM)+PINACT) GO TO 820 . 

ACTIVE CONSTRAINTS 
TYACT is no. OF TYPES ACTIVE 
KACT IS POINTER TO ACTIVES 
K=K + 1 

TYACT-(K)=TYPECI,KSYM) 

KACT(K,KSYMD=KPTSCI,KSYM) 

SMLCK',KSYM]=ST8M(I,KSYM) ' 

ACTIVE Cl fKSYM) = ACTDES(l) 

GO TO 830 

INACTIVE CONSTRAINTS 
nitcksymi=nit-cksym) + i 

NQsNITCKSYMD 

KINACTCN0,KSYM)=KPTS(I,KSYM)- 

ACTlvECIfKSYM)=ACT0ES(2) 

CONTINUE 

IF ( (KEYOUT(O).EO.n .AND. ( CkPR . EQ. KPR INT ) . OR . ( I TER . EQ . KST ART J . OR . 

1 ITER.EQ.KQUIT)) ) CALL OUTPT (2,KSYM,KP0INT/CT C1,KSYM) ,OMEGA, ITER, 
2,KPTS(lf KSYM) ,KSM,STBMC1 ,K3YM) ,RQ(1,KSYM) ,TYPE(1 .-KSYM] ,ACTIVE(1/K; 
3YM).,KIN,KOUT,2A,Z8,ZC/ZD/ZE,PA,PB,PC,PD,PE/Ni,N2,MWM2,GAIN,KONT, 
«,BrC,KOLD) 

NML (KS.YM)=K 

CALL PARTAL (OMEGA, NML CKSYM} ,CT Cl/ KSYM) , KACT n ,KSYf-n ,TYACT(1) ,T,G 
1GC,KSYM,C1 ,CI,W0RKI,A,B,C,D,£,F,C5,PG(NAM+1 ,n ,PGCNAM.+ 1 ,KN0T + 1),K 
2N,K0UT,ZA,ZB,ZC,ZD,ZE,PA,!^S,PC,PD,PE,N1 ,.N2,M1,M2,GAIN,K0NT) 
NAM=NAM+NML(KSYM) 

CONTINUE 

normalizing' PG 

DO 860' 1 = 1 ,NAM 
SUM=0. 

DO 850 J=1,NPARC 
•SUM=SUM+PG(I, J3**2 . 

SUM=DSQRT(SUM5 
DO 860 J=1,NPARC 
IF CSUM.LT^O. lE-06) SUM=1.E*06 
PGCI,J0=P6(I, JO/SUM 
IF CITER.GE.KOUIT) WRITE (6,870) 

FORMAT C'0',5X,'*** TERMINATION REASON: MAXIMUM ITERATIONS ***') 

IF (ITER.GE.KGUIT) GO TO 900 
IF (NAMI.NE.O) go to 920 

write (6,880) 

FORMAT C '0 ', 15X, '**** ALL SYSTEM REQUIREMENTS HAVE BEEN MET ****') 
WRITE (6,890) ITER 

format ('0*,5X,'THE SUBOPTIMAL compensator OCCURRED ON ITERATION N 
10.; ',I3,/5X,'W1TH THE FOLLOWING COEFFICIENTS; '/) 

DO 910 KSYM=i,KIN 
N=NM(KSYM) 

IF (KEYOUTU ) .EQ. ll CALL OUTPT ( 3 , KS YM , KPO I NT , CT ( 1 , KS YM) , OMEGA , I TE 
1R,N,KPTS(1,KSYM),KMINCKSYM) ,STBM(1 ,KSYM) ,RQ(1 ,KSYM) ,TYPECi ,KSYM) ,A 
2CTIVE(1,KSYM),KIN,K0UT,ZA,ZB,ZC,ZD,ZE,PA,P8,PC,P0,PE,N1 ,N2,M1 ,M2,G 
3AIN,K0NT, A,B,C,KCLD) 
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STOP 

SET DOT PPOOUCT 
DO 930 N=1,NAM 
WElGHT(N)sl, 

IF (NAM.LE.LPV) GO TO 960 
WRITE (6,950) 

FORMAT ('O', 5X, '***** TERMINATION REASON - NO. OF ACTIVE CONSTRAJ 
INTS EXCEEDS THE NO. OF ALL0^>A8LE VARIRA8LES *****') 

STOP 

CALL DIRVEC (PG, NAM, NPARC,DV, WEIGHT, PFX, PFY,W0RKR,E,F,C5,E2) 

KRE = 0 
KKK = 0 

CALL XCHECK (K IN, KOUT , ZA , ZB , ZC , ZD, ZE , N1 , N2 , D V , ZET AZ, PG , NAM, KRE, NZE 
'1R01,NZER02,LPV,A,B,C,0,E,F,KKK) 

IF C(KR£.EQ.3) .ANO, (ITER.EQ.KSTART) ) GO TO 1000 
IF (KRE.EO.l) GO TO 9-30 

CALL XCHECK (KIN,K0UT,PA,PB,PC,PD,PE,M1,M2,DV,ZETAP,PG,NAM,KRE,NPC 
1LEI,NP0LE2,LPV, A,8,C,D,e,F,KKK) 

IF C CKR£,£Q.3).AND, (ITtR.EQ.KSTART) ) GO TO 1000 

IF (KRE.EQ.l) GO TO '9A0 

PSQ=0, 

DO 970 Isl,NPARC 
PS0=PSQ+DVCI)**2 
PMG=S0RT(PSG1) 

IF (PMG.LT.l.E-08) PMG=l.E-08 
DEL=STEP/PM6 

KKK = 1 

CALL CHANGE ( Z A , ZB , ZC , ZO, ZE , N1 , N2 , 0 V , DEL , KKK , K IN, KOUT , A ,'8, C ) 
call CHANGE (PA,PB,PC,P0,PE,M1 ,M2,DV,DEL,KKK,KIN,K0UT,A,8,C) 

IF (KRE.NE.3) STPOLD=STEP 
IF (KRE.e0.3) STEPsSTEP+STPOLD 
KKK = 0 

IF CKRE.E0.3) GO TO 990 
KRE = 0 

CAUL YCHECK (K IN , KOUT, M ,N2, ZA , ZB , ZC, ZO , ZE, KRE, STPOLO, PMG, ZET AZ , NZ 
lEROl ,NZER02,DV,A,6,C,D,.E,KKK,STEP) 

CALL YCHECK (KIN, KOUT, M1,.M2,P A, PB, PC, PD, PE, KRE, STPOLD, PMG, ZETAP,NP 
10LE1,NP0LE2,DV,A,B,C,D,E,KKK,STEP) 

IF (KRE.EQ.3) GO TO 980 
KRE=0 

ITER=ITER+1 
GO TO 370 
WRITE (6,1010) 

FORMAT ( '0',5X, *•*** TERMINATION REASON: INITIAL COMPENSATORS DO N| 

lOT SATISFY ZETA CONSTRAINTS ***') 

STOP 
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SURROUTIME ADOPTS 

SUBPOUTINE adopts ( KPC) INT » kIM/ KOUT , rjS , STRM, kPTS / CT , G/ GC / T , OMEG 
1 ,KG08Ak,KPTWAX,NIT,KIA!ACT/NML,KACT,KP0LP/KC1LD, KSYN^Cl ,CI,»vOPKI , ZA 
2ZB, ZC,ZD,ZE,PA,P3,PC,PD,PE,M1 ,N2,^''1 ,M2,GAIt^),K0NT, A,e,C,D,E) 

SUBPROGRAM DESIGNED TO GENERATE ADDITIONAL FREQUENCY DATA; 

PROGRAM INCORPORATES THE ROUTINES CRT, INTER, AND TRFR. 

SUBPROGRAM VARIABLES! 

KPOINT « INTEGER NUMBER OF CURRENT DATA POINTS 

KIN - INTEGER NUMBER OF CONTROLLER INPUTS 

KOUT - INTEGER NUMBER OF CONTROLLER OUTPUTS 

N8 - starting NO, OF MARGINS TO BE INVESTIGATED 

NM - INTEGER NUMBER OF MARGINS INVESTIGATED 

STBMd) - REAL ARRAY OF STABILITY MARGINS 

KPTS(I) - INTEGER array OF FREQUENCY NOS. OF MARGINS 

CTCII - COMPLEX ARRAY OF FREQUENCY RESPONSE 

G(I,J,K5. - 3D COMPLEX ARRAY OF DISCRETE FREQ, RESPONSE 

GCCI,J,K)- COMPLEX ARRAY OF COMPENSATION EVALUATED AT DISCRETE PT 

T(I,J,K) - complex array of TRANSFER RESPONSE Gi.*G IS STORED 

OMEGA (I) - complex ARRAY OF DISCRETE FREQUENCY POINTS 

KPTMAX - MAXIMUM NUMBER OF FREQUENCY POINTS ALLOWABLE 

NIT ” AN INTEGER OF INACTIVE MARGINS 

KINACTCD- ARRAY OF INTEGERS CORRESPONDING TO INACTIVE MARGINS 
NML - AN INTEGER DENOTING ACTIVE MARGINS DETECTED 

KACT(I) - FREQUENCY DATA NUMBERS OF ACTIVE MARGINS 
KPOLD - INTEGER OF DATA POINTS OF LAST ITERATION 

KOLOCI) - INTEGER ARRAY OF PREVIOUS DATA POINTS 

KSYM - INTEGER REFERENCE TO PARTICULAR SUBSYSTEM 

COMMON ITER,KSTART 
INTEGER A,8,C,D,E 

complex ct,g,gc, omega,? 

DIMENSION CKO, A}, G(8,A,D), GC(A,8,D), KACKE,A), KINaCT(E,A), K 
ILDCn, KP7S(E,A), NIT(A), NM(A), NML(A), QMEGA{0), STaMCl), T(A,A, 
20 ) 


KGOBAK=0 

NX=NM(KSYM) 

DO no N=N8,NX 
K=KPTS(M,KSYM) 

OG 1=CABS(CT(K, KSYM) -C T ( K“1 , KSYM) ) 
DG2=CABS(CTCK+1 ,KSYM)-CT (K,KSYM) ) 
IF (ST8M(M) ,GT,5,*0G1 ) GO TO 100 
K=K + 1 

IF (K,£Q,2) GO TO 100 

CONTINUE 

KPOINT=KPOINT+l 

IF (KPOINT, GT. KPTMAX) RETURN 
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SEPARATING MATRICES TO ADD NEW FREQUENCY POINT 
00 30 M=KPOINT,K,-l 
OMEGA{M)=OMEGA(M-1) 

‘ 00 20 1=1, KIN 
CT(M,I)=CT(M-1,I) 

DO 20 J=1,KIN 
TCI, J,M)=T(1, 

DO 30 1=1, KIN 
DO 30 J=1,K0UT 
G(J,I,M)=G(J,I,M-n 
GCCl,J,M)=GC(I,J,M-n 
KG08AK=1 

CALL INTER (OMEGA (K-1 ), OMEGA (K-2) , OMEGA (K-ll ) 

00 «0 I=1,K0UT 
DO «0 J=1,KIN 

CALL INTER (G C I , J, K- 1 ) , G (I , J ,K-2 1 , G ( I , J , K-t) ) 

CALL EVAL (OMEG A (K-1 ) , GC , G, T , K-1 , K IN , KOUT , Z A , ZB , ZC , 2D, ZE, PA, PB , PC , 
1PD,PE,N1 ,M2,M1,M2,GAIN,K0NT, A,B,C,D) 

00 50 I=1,KSYM 

CALL CRT a,K-l,I,TU,l,K-l),KIN,KOUT,CT(l,n,G,PCG,Il,Jl,Cl,CI,W 
1RKI,A,8,0) 

NX=NMCn 
00 50 L=1,NX 

IF (KPTS(L,I).GE.K-n KPTS (L, II =KPTS (L , 1 1 + 1 
IF (ITER.EQ.KSTART) GO TO 80 
DO 70 I-l/KIN 
NYsNITd) 

NZ=NML(I) 

DO 60 L=1,NY ■ 

IF (KINACTCL,I).GE.K-n KINACKL , 1) =KINACT (L, I ) +1 
DO 70 L=1,NZ 

IF (KACTCL, n ,GE.K-n K ACT (L , IJ =KACT (L, I ) +1 

CONTINUE 

DO 90 L=l,KPOLD 

IF (KOLD(LKGE,K-n K0LD(L)=K0LD(L) + 1 
IF (ST0M(N) ,GT.5.*0G2) GO TO 110 
STBM(N)=6.*DG2 
K=K+2 

IF (K.GT.KPOINTl GO TO 110 

GO TO 10 

CONTINUE 

RETURN 
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SUBROUTINE CHANGE 


’SUBROUTINE CH ANGE (Z A » ZB , ZC / ZD , ZE> N1 , N2, OV» DEL » KKK , )< IN, KOUT , A, 8f C ) 

|C c . 

Ic OE-rrIGNEO TO CHANGE COMPENSATOR COEFFICIENTS ACCORDING TO THE DIRECT 
'C lONAL VECTOR 
C 

INTEGER A,8,C 

DIMENSION ZACA,8,C),28(A,B,C),ZCtA,B,C)>ZD(A,e,C),ZE(A,B/CJ ,OVU) 
C ,M (A,B) ,n2(A,BI 

on <1 islfKIN 
DO ii Jsl.KOUT 
NV = 'Jl(I,J) 

IF(NX.EQ.O)GO TO 2 
DO 1 L*l,NX 

2A(I,J,L)=ZACI, J,L)+DV£KKK)*D£L 
ZBCr, J,L)=Z“ n, J,L)*DVCKKK + n*DEL 

1 KKi< = KKK+2 

2 NXr'M2(I,J) 

IF(NX.EQ.O) GO TO « 

00 5 L=1,JX 

2Cn,J,U=ZC(I, J,L)«’DVCkKK)*D£L 
20(1 , Jr Ulszon ,J',L)+r.VCKKK + n*OEU 
ZE(I, J,U=ZEn, J,L)+0 V(kkK + 2)*DEL 
S KKK = t<Krt + 3 

tt continue 

return 

END 


SURROUTINE CRT 


SUBROUTINE CRT (RE Y r K , KSTM, T, KIN, KDUT,CTr P, PCG, 11 , J I ; C I , C I , v«.ORKI , 

subprogram performs 2 FUNCTIONS DENOTED BY THE INPUT KEY: 

KEY IS I, CLOSED LOOP FREOUENCV RESPONSE CT(S) IS FOUND IN 
TERMS OF INPUT VECTOR; KEY IS 0, SUBPROGRAM AIDS PARTAL IN 
DETERMINirjG THE PA'RTIAL OF CCS) wRT. G(S)IJ. PROGRAM IN- 
CORPORATES MATTTiC, matmuL ROUTINES, 

SUBPROGRAM VARIABLES; 

K£V - PROGRAM MODE VARIABLE 

KSTM - SUBSYSTEM IN COf.SIDEPATION 

T(I,J) - complex transfer RESPONSE G(S)*P{S) 

KTJ - NO, OF CONTROLLER INPUTS 

KOUT - NO. OF CONTROL OUTPUTS 

P(I,J) - COMPLEX PLANT RESPONSE ARRAY 

II - CONTROL INPUT INDEX OF G(S)IJ 

J1 - CONTROL OUTPUT INDEX OF GCS)IJ 

C T ( I) - C QMPLEX CLOSED LOOP FREQUENCY RESPONSE 



c 

C TK IS -THE KTH COLUMN OF T(S) 

C Cl IS THE KSTM RESPONSE 

C ZEROING THE KTH COLUMN OF (T(S) + I): 

C 

INTEGER A,B,D 

complex C1(A,A),T(A,A),CI(A,A),P(B,A),CT{0)»PCG 
DO 10 I=1»KIN 
DO 10 J=l,KIN 
Cl (I,JI=GMPLX(0.,0.) 

IF CL.NE.KSTMI C1CI,J)=T(I,J) 

IF (I.EQ.KSTMD C1(I,J)=CMPLX'C0,,0.) 

IF (I.EG.J) Cl Cl, JDsCwPLXCl ,,0.)+Cl (I»J) 

10 CONTINUE 

call MATINC (C1,KIN,CI,I£R,IvORKI,A,A,2*A) 

IF (IER.E0.2) GO TO «0 

call MATMUL (T,CI,C1,CI,T,T,KIN,KIN,KIN,KIN,I,0,A,A,U 
CT(K)=C1 (KSTMrKSTM) 

C Cl AT THIS POINT IS total RESPONSE C IN NOTES. 

IF (KEY.EQ.l) RETURN 
C NOl^ THE PARTIAL OF C W.R.T, G(I,J) 

PCG = 0-.0 
DO 20 1=1, KIN 

20 PCG=PCG+P(J1 ,I)*CI (T.KSTM) 

IF (11 .EO.KSTM) GO TO 30 
PCG=-C1(KSTM,II)*PCG 
30 CONTINUE 

RETURN 

aO aPITE C6,50) KSTK,K 

50 FORMAT (/bX, 'INVERSE MATRIX INDETERMINATE BY THIS METHOD AT KSYM= 

l',I3,'K= ',13/) 

RETURN 

C 

END 


C SUBROUTINE DELETE 


SUBROUTINE OELETE(KPOINT,kIN,KOUT,OMEGA,G,NIT,KINACT,NML,KACT, 
CIJER,KPOLO,KOLD,KNEh,A,B,C,0,E) 

INTEGER A,B,C,D,E 

complex omega, G 

DIMENSION 0MEGA(D)',G(B, A,D),KINACT(E,.A),NIT(A) ,KACTCE, A) ,KOLD(D), 
CNMLCA),KNEW(D) 


L = 0 

DO 1 K=1 ,KPOINT 
DO a 1=1, KIN 
NX=NIT(I) 

IFCNX.EQ.O) GO TO 2 
DO 1 J=1,MX 

IFCK.EQ.KINACTCJ,!)) GO TO 6 
NXsNML C I ) 

IF(NX.EQ.O) GO TO U' 


il 

2 
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DO S J=1,NX 

3 IF(K.EB.KACTCJ,n) 60 TO 6 

:a CONTINUE 

DO s J=1,KP0L0 

5 IF(K.EQ.KOLDCJ)) go TO 6 
GO TO 7 

6 L=L+1 
KNEW(U=K 

7 CONTINUE 
aD£LET=KPOINT-L 
KP01NT=L 

do' 8 K=1,KP0INT 
L=KNE'.v(K) 
riMEGACK) = OMEGA(L) 

00 8 I=1,K0UT 

DO 8 J=1,KIN 

8 G(I.J»i<)=G(I,J/L) 

DP 12 K=1,KP0INT 

00 12 1=1, KIN 

NXsNITCI) 

00 9 J=1,NX 

9 IF(K-vE>v(K),EQ.KlNACTCJ,inM'NACT(J,n=K 
riXsMNLCn 

DO 10 J=UNX 

10 IFCKNEi';(K) .EQ.KACT(J, 1) )KACTCJ,I)=K 

no 11 J=1,KP0LD 

11 IF (K0LD( J) .E‘3.KNEw(KnK0LD(J)=K 

12 CONTINUE 

•V9ITE (fe, 13)K0ELET,ITER 

15 FDP'1ATC1H0,5X, 13, 2X, 'POINTS ivERE DELETED ON ITERATION NO, '13) 

RETURN - 
C 

end 


C SUBROUTINE OIRVEC 


subroutine DIRV£C(G,NN,KPARC,0V,WEIGHT, A1, ai,workr,e,f,z,h) 
directional vector program 

subprogram designed TO- CALCULATE THE DIRECTIONAL VECTOR OF 
the constraint IMPROVEMENT ALGORITHMM. DIRECTIONAL VECTOR 
DV CALCULATED AS 

DV = («G)*A 

WHERE )“G IS -CN,M] GRADIENT MATRIX WHOSE COLUMNS ARE THE 
GRADIENTS OF THE ACTIVE CONSTRAINTS; COLUMN VECTOR. A IS 
- ‘ A = INVCi!>G'«(G)*C; 

.COLUMN 'C rs M- COMPONENT VECTOR OF POSITIVE ELEMENTS; M-'lS 
NUMBER OF COMPENSATOR COEFFICIENTS, 


I* 
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DEFINITIONS OF I/O VARIABLES 

G -matrix rtHOSE ROWS CONTAIN THE GRADIENT VECTORS OF THOSE 
STABILITY margins ONLY CONSIDERED PERTINENT 
NM -number of stability margins CONSIDERED PERTINENT 
KPARC -NO. ROWS IN G, I.E., COMPENSATOR COEFFICIENTS 
DV(I) -real array of directional VECTOR 
WEIGHT-WEIGHTING FACTOR VECTOR 

integer E,F,Z,H 

DIMENSION G(E,F),,A1 CE , E ) / A I (E» E 3 , WE IGHT ( 1 ) , OV ( 1) 

DO 2 K=1,NM 
DO 2 J=K,NM 
SUM=0. 

DO 1 1=1, KPARC 

1 SUM=SlJM+G(J,I)*G(K,n 
Al CJ,K3= SUM 

A1 (K, J)= SUM 

2 CONTINUE 

IFCNM.GT.n GO TO 3 
Ain, 13 = 1. /Aici,n 
Al (1,13= v-EIGHT{n * AIC1,13 
GC TO b 

3 CONTINUE 

call M a TJ N V.( Al , NM , a I , I E R , WORK P-, E , E , H 3 
I'f(IER.EQ.03 go TO .5 
rtRITE(6,«3. , . . ^ 

A FOR'-'AK'O*’, lbX„ '**** TERMINATION: PAR.TI'ALS NOT LINEARLY INDEPENO 

CNT ****'3 
STOP 

5 CALL MATKUL(AI,Y,Y,AI,w£IGHT,A1,Nm,nm,NM,1,1,1,E,E,13 

b CONTINUE 

•00 3 1 = 1, KPARC 
SU^/s ,0. 

DO 7 .J=1,NM 

7 SUM=SUM-*-G(J,n*Al (J,13 

8 DV(I3=SUM 
RETURN - 

END 


C SUBROUTINE EVAL 


SUBROUTINE EVAL ( XF » GC , Gf T , K, KIN, KOUT , ZA , ZB , ZC , ZO, ZE , PA , PB , PC, PD, 
lfc,Nl,N2,Ml ,M2,6ArN,K0NT, A,B,C,D3 

INTEGER A, 8', C/D 

complex GCX,GCY,GCN,GCD,GC(.A,B,D3 ,G(B,A,DD ,TCA,A,D3 
DIMENSION N1(A,B3, C0EF(33, ZA(A,B,C3, ZBCA,B,C3, ZC(A,B^C3, ZD(A 
1B,C3, ZE{A,e,C3, PA(A,R,C3/ PB(A,B,C3, PC(A,6,C3, PDCA,B,C3, PE(A 
2B,C), GAIN(A,B3, N2(A,B3, MHA,B3, M2(A,83 
COMPLEX XF 




c EVALUATION OF GC r THE COMPENSATOR TRANSFER 

DO '>0 1 = 1, KIN 
DO 90 J=l,KOUT 
GCX=CMPLX(1 , ,0.) 

GC>=CMPLX(1. ,0) 

GCN=CMPLX(1.,0.) 

GCD=CMPLX(U ,0.) 

NCsNl Cl, J) 

IF (NC.EQ.O) GO TO 20 
DO 10 L=1,NC 
C0EFC1)=ZA(I,J,L) 

COEF (2)=ZB(I,J,L) 

CALL POLEV (COEF,l,XF,GCX) 

GCN=GCX*GCN 
10 CONTINUE 

20 MC=M1(I,J) 

IF C*^C.EQ,0) GO TO aO 
DO 30 L=1,MC 
COEFCl )=PACI,J,L3 
CriEF(2)=PBCI,J,L) 

CALL POLEV (COEF, 1 ,XF,GCY) 

GCD=GCV*GCD 
30 CONTINUE 

«0 IF CGCD.EO.O.] GC0=CMPLX(1.,0.) 

GCN=GCN/GCD 
'4C = N2(I,J) 

IF (NC.EO.O) GO TO 60 
DO 50 L=1,NC 
CnEF(,n=ZC(I ,J,L) 

COEF{2)=ZO(I,J,L) 

COEF(3)=ZE(I,J,L) . 

call POLEV CCOEF,2,XF,GCX) 

GCN=GCX*GCN 
50 CONTINUE 

60 mc=m2(1,J) 

GCDsCMPLXCl. ,0.) 

IF (MC.EQ.03 60 TO 60 
DO 70 L=1,MC 
COEFCn=PCCI,J,L) 

C0EF(2)=P0(I,J,L) 

C0EF(3)=PE(I, J,L) 

call POLEV (C0EF,2,XF,GCY) 

GCD=GCY*GCO 
70 CONTINUE 

80 CONTINUE 

IF (GCD.EO.O.) 6CD=CMPLXC1.,0.) 

GC Cl, J,K1=GAIN(I,J)*6CN/GCD 
90 CONTINUE 

C GC'*G transfer FUNCTION 

C HERE STARTS TRANSFER GC*G O.L. FREQUENCY RESPONSE 

call MATMUL CGC(l,l,K),GCl,l,K3,T,GC,GC,GC,KIN,KOUT,KOUT,KIN,K,0,i 

1 , 8 , 8 ) 

RETURN 

C 

end 
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C SUBROUTINE GAINMG 


SUBROUTINE GAINMG (GTOTAL# KPOINT, NM, FOMIN, FQMAX , KPTS/ STBM, OMEGA) 
C SUBPROGRAM FOR CALCULATING GAIN MARGINS 


C definitions of I/O VARIABLES 
C 

ic GTOTAL-COMPLEX array OF COMPENSATED OPEN FREQ. RESPONSE 
jC KPOINT-NO. of POINTS 
jc OMEGA -ARRAY OF FREQS, 

C NM -COUNTER 

C KPTS -FREQUENCY NOS. WHERE .MARGINS' OCCUR 
p ST6M -stability, margins OF MARGINS 

C FQMIN -LONER FRE'O. FOR KARGIN DETECTION 

C FQMAX - UPPER FREQ, for MARGIN DETECTION 

complex 0'*EGA,GTDTAL 

, DIMENSION GTOTALU), KPTS(l), STBM(l'), OMEGACl) 

Psl.O 

DO '50 I = 1,KP0INT 
SBsAIMAGCGTOTALCD) ■ 

IF CI.E.Q.1) S1=S2 

■ IF {AIMAGtOMEGACn ) .GT.FQKAX) RETURN 
IF C'AIMAGCOMEGACm.LT. FOMIN) GO TO 20 
IF (ARS{S2) .LT.1 .OE-20) GO TO tO 
SGN=S2/ABS(S2) 

IF (S1*SGN.GT.0.0) GO TO 20 

10 IF C(realcgtotalci)).ge.o.).and.(real(Gtotal(i«i)).ge.o.)) go to 

10 

Ii=I-l 

IF CA8S(S2),LT,ABSCS1)) 11=1 

NM=NM+1 

KPTS(NK)=I1 

FKAC=S1/(S1-S2) 

STBMCNMI=CABSCP+FRAC*GTOTALn )+(l.-FRAC)*GTOTAL(I-l)) 

20 Sl=S2 

30 continue 

RETURN 


END 
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SUBROUTINE INTER 


subroutine INTER(S,T,R) 

INTERPOLATION SUBPROGRAM 

SUBPROGRAM INTERPOLATES SPECIFIED INPUT DATA IN 
CONJUCTION WITH ADOPTS ROUTINE; MAGNITUDES ARE 
interpolated logarithmically, PHASES LINEARLY, 

subprogram VARIABLES; 

S - COmPLEY L0i’>ER BOUND OF QUANTITY ‘ INTERPOLATED 
T - COMPLEX UPPER BOUND OF QUANTITY INTERPOLATED 
R - complex resultant of THE INTERPOLATION 

COMPLEX S,T,R 
DATA PI/3. 1«1BR26535697R/ 

C 

X = SQRT(CABSCS)*CABSCn) 


Y = 0. 

IF((x.GI.l.E-lOn GO TO I 
GO TO 2 

1 CONTINUE 

U = ATA\12(AIMAG(SI #REAL(Sn 
V=ATAN2(AIMAGCT) ,REAL(T)) 

IFCU. LT.O.) V=V+2.*PI 

IFCV. LT.O.) V=Vf2,*PI 
Y=CU+V)/2. 

IF C ABSCU-Vl.GT.PI ) Y=Y-PI 

IF(AHSCY-U).GT.PI/2.)Y=Y-PI 

2 R=CMPLX(X*COSCY) ,X*SIN(Y)) 
RETURN 

C 

END 




SUBROUTINE MATING 


SUBROUTINE MATING (XX,N, YY, IER,X»A,B,.G) 

SUBROUTINE FOR FINDING AN INVERSE OF A MATRIX 

X = MATRIX FOR WHICH INVERSE IS .TO BE TAKEN 
^4 = dimension of square MATRIX X 
Y = INVERSE OF X 
lER - ERROR CODE 

IEH= 0 - NO ERROR EXISTS 

IER= 2 - MATRIX DOES NOT POSSESS AN INVERSE 
INTEGER A,8,C 

DIMENSION XX(A,4), YY(A,AJf X(A,C) 

I'^PLICIT complex (D-H,0-Z) 


DO 10 1=1, N 
DO 10 J=1,N 
X(I,Jj=XX(I,J) 

|10 CONTINUE 
IER=0 
N2=N+N 
N1=N+1 
00 20 1=1, N 
DO 20 K=M,N2 
L = K-N 

IF (L.NE.1) X(I,K)=CMFLX(0.,0.) 

IF (L.eo.n x(i,k3=cmplx(i.,o.) 

20 CONTINUE 

PC bO 1=1, N 

IF CCA8S(XtI,n).LT.1.0E-25) GO TO 70 
SCRv\=X(l , I) 

DO 30 K=1 ,N2 

30 XU ,K) = y(I,K)/SCRrt 

00 50 L=1,N 
IF (L.EQ.U GO TO 50 
SLOP=X(L,I) 

DO <40 K = 1,'N2 

XCL,K)=X(L,K)-SLOP*X(I,K) 

UO CONTINUE 

50 CONTINUE 

DO 60 1=1, N 
00 60 K=N1,N2 
L=K-N 

60 YY(I,L)=X(I,K) 

GO TO 6b 
70 IER=2 

60 RETURN 

C 

END 
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SUBROUTINE MATINV 


subroutine MATINV ( XX , N , YY f lER » X,E , Zr H) 
INTEGER E^Z#H 

DIMENSION XXCE^Z), YY(E,Z), X(EfH) ' 

DO 10 1 = 1, N 
DO 10 J=1,N 
10 XCI,J)=XX'(I,J) 

N2=N+N 
NlcN+l 
DO 20 1=1, N 
DO 20 K=N1,N2 
L=K-N 

IF CL.NE.I) XCI,K)=0. 

IF CL.EQ.n X(I,K) = t, 

20 CONTINUE 

DO SO 1 = 1, N 

IF (ABS(XCI,n).LT.l,E-25) GO TO 70 

SCRvv=X(I,n 

DO 30 N=1,N2 

30 X(I,K)=X(I,KI/SCftW 

DO SO L=1,N 
IF (L.EQ.U GO TO 50 
SLOP = X(L,n • 

DO UO K = 1',N2 

X(L,K)=X(L,K)-SLOP*XU,K) 
ao CONTINUE 

50 CONTINUE 

DO 60 1 = 1, N 
DO 60 K=N1,N2 
L = K-N 

60 YY(I,L)=X(I,K) 

GO TO 80 
70 IER=2 

80 RETURN' 

C 

END 


C SUBROUTINE MATMUL 


SUBROUTINE MATMUL C AC, BC , CC , AR , 3R, CR , N , L » LI » M, NO , NC , A,B, C ) 

SUBROUTINE MATRIX MULTIPLICATION 

MATRIX MULTIPLICATION MULTIPLIES A * B 
MATRIX AC IS N X L 
MATRIX 8C IS L X M 

matrix CC IS The RESULTANT MATRIX N X M 
MATRIX CR IS REAL RESULTANT 
NC IS complex KEY; 

NC=0; AC,BC,CC, ARE COMPLEX MATRICES 
NC=1F AR,8R,CR,ARE REAL MATRICES 




INTEGER A,fl,C 
CO»-'PLex AC,BC,CC 

DIMENSION AC(A,B), BCCB,A), CC(A,A,NO), ARCA/B), BR(B,C), CR(A,Cf 
10 ) 

C 

IF CNC.EG.13 GO TO 20 
DO JO 1=1, N 
DO 10 J=1,M 

CC(I, J,ND)=CMPLX(0.,0.) 

DO 10 K=1,L 

10 CC{I,J,ND)=CC(I,J,ND)+ACn,K)*BC(K,J) 

return 

20 DO 30 1 = 1, N 

DO 30 J=1,M 
CR(I,J/ND)=0, 

DO 30 K=1,L 

50 CRCI; J»N0)=CRCr, J,ND)+AR(1,K)*BRCK,J) 

RETURN 

END 


SUBROUTINE NYOIST 


subroutine nyQIST CN,G,NRHP,NC0N,NCIRL,MZ,FMAX,F) 

COMPLEX GCn,F(n 

DOUBLE PRECISION PI 

data pi /3. 1J159265358979D+00/ 

KC=tJCON/2 

LC=(NCON+l)/2 

Cl=ATAN2CAIMAGCGm)., l.O+REALCGU))) 
SUM=-CUA8S(ATAN2(0.0,1.0 + REAL(GC1)) ))*ABS(Cn/Cl 
IF CKC.NE.LC) SUM=0.0 
DO 10 1=2, N 

C2 = ATAN2{A1MAGCG(I)), l.OtREALCG(I) )) 

0IFE=C2-C1 

IF (ABS(DIFF) .GT.PI) DIFF=C2-C1+2,0*PI*ABS(C1)/C1 
SUM=SUM-OIFF 

IF (AIMAGCFd) ) .GT.FMAX) GO TO 20 
0 C1=C2 

I=N 

0 SUM=+C2-ASS(ATAN2(0.0,1.0+kEAL(G(I))) )*ABS(C2)/C2+SUM 

SUM=2.0*CSUM+C2)+P1*NC0N 
SUM=SUM+ABS(SUM)*PI/(ii.O*SUM) 

NCIRL=SUW/(2.0*PI) 

NZ=nRHP+NCIRL 

RETURN 


END 
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SUBROUTINE OUTPT 


subroutine OUTPT (N r KSVM, KPQINT » CT , OMEGA , ITER f i^M/KPTS , KMIN, STBM» RQ 
1 / TYPE, ACTIVE,KIN,KOUT,ZA,ZBf ZC,ZD,ZE,PA,PB,PC,PO,PE»Nl ,N2,M1 ,M2,GA 
2IN,K0NT, A,8,C/K0L0) 

SUBPROGRAM DESIGNED TO OUTPUT INFORMATION IN THREE AREAS 

AS designated by the key N; 
key: 

N=0, OUTPUT compensator INFORMATION 
N=l, OUTPUT FREQUENCY RESPONSE 
N=2, OUTPUT, STABILITY MARGINS INFORMATION 
N=3, complete OUTPUT INFORMATION 

INTEGER A,B,C 

INTEGFR type , ACTIVE, XXP, ZOUM 0 
COMPLEX CTCn ,DmeGA(1),)( 

dimension GAIN(a,B), KDNTCA.B), NJ(A,B), N2(A,8I, M1(A,B), M2(A,B) 

1, ZA(A,B,CI, ZB(A,B,CI, ZCCA,P,C), ZD(A,8,C), ZE(A,6,C), PA(A,8,CI 

2, PB(A,B,C), PCCA, B,C3, PD(A,B,C), PECA,8,C), KPTS(l), STBMd), RQ 
3(13, ACTI-VECl), TYPEd), KOLDCl) 

DATA lbK,IAT,XXP,RA62 /IH , 1 H* , I HP , U « . 59 1 559/ 

IF {N.EQ.23 GO TO 170 
IF (N.EQ.l) GO TO 110 
IF (KSYm.eq,!) i.,ritE (6,20) 

DO 100 1=1, KIN 
DO 100 J=1,K0UT 

IF (gain(i,j).eo.o.) go to 80 

iVRITE (6,30) I, J.,GA1N(I,J) 

NCsNl (I,J) 

IF (NC.NE.O) V.RITE (6,90) ( Z A ( I , J , L) , ZB (I , J, L ) , L = 1 , NC ) 

NC=N2d,J) 

IF (NC.NE.O) WRITE (6,50) t ZC ( I , J , L) , ZD ( I , J, L) , ZE ( I , J, L) , L = 1 , NC) 

MC = Mld,J) 

IF (MC.NE.03 WRITE (6,60) ( P A ( I , J , L) , PB ( I , J , L) , L= 1 , MC ) 

MC=M2(I,J) 

IF (MC.NE.O) WRITE (6,70) (PC C I , J, L) , PD ( I , J , L ) , PE ( I , J , L) , L=1 , MC) 
WRITE (6,10) I,J,KONT(I, J) 

FORMAT ('O'fSX,' DC GAIN CONSTRAINT FOR EACH CHANNEL ( IF KONT=l , 

1 ALLOWED TO VARY; -IF KONT =2 , HELD CONSTANT )”,/, 7X , 'KONT ( ',13,', 
2'. 13,') = ',13,/) 

FORMAT {/,5X,'THE COMPENSATION ELEMENTS ARE DESCRIBED BY TRANSFER 
1 FUNCTIONS IN', /,5X, 'CASCADED FIRST AND SECOND ORDER FACTORS'; ',/, 3 
21X, 'M ',mx, 'n2',/31X, 'PR -(2A + ZB S) PR (ZC + 2D S + ZE S**2)', 

3/,30X, '1 = 1 ',aX, 'I',5X, 'I',3X, 'J=l ',3X, 'J'5X, 'J'6X, 'J',/18X, '(GAIN) 

^ — - .',/,31X,'M 

51 ' , 19X, 'M2',/, 31X, 'PR (PA + PB S) PR (PC + PD S + PE S**2)',/,30 
6X,'I=1 I 1 J=1 J J J'/) 

FORMAT C/,5X,'COMPFNSATOR(',I3,'r',I3,'3; GAIN = ',G10.5,/) 

FORMAT (15X, 'COMPENSATOR COEFF IC lENTS; ' , / , ( 1 5X , 'Z A = ',G12.6, lOX, ' 
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*COMPEMSATOR(',I3»'/'rI3,') HAS ZERO CONTRIBUTION.') 


IZB = ',G12.6)) 

FORMAT (15X,'ZC s ',G12.6, lOX, 'ZD = ',G12.6, lOX, 'ZE s ',Gl2.6) 
FORMAT C15X,'PA = ' , G 1 2 . 6, 1 OX , 'PB s '»G12.6) 

format (15X,'PC = ',G12.6, lOX, 'PD = ' , G1 2. 6, 1 OX , 'PE * ',G12,6) 

GO TO 100 
hRITE (6,90) I,J 

FORMAT (lOX, 'COMPEMSATOR(',I3,',',I3,') HAS ZERO CONTRIBUTION.') 
CONTINUE 

IF (N.NE.3) RETURN 
hRITE (6,120) KSYM 

FORMAT (/,5X,'THE FREQUENCY RESPONSE OF SYSTEM NO. ',13,' IS: ',/ 
17x, 'DATA', lOX, 'COMPLEX OMEGA ', 1 0 X , 'COMPLEX CT',10X,'MAG & PHASE'# 
2 ) 

L=I 

DO laO Krl,KPOINT 
ZDlJMBsIBK 

IF (K.NE.koLO(L) ) GO TO 130 
ZDUMHsI AT 
L = L + 1 
CONTINUE 

X=CMPLX (CABS(CT(K) ) , ATAN2(AIMAG(CT(K) ), REAL (CT (K ))) *57 .295779) 
tvPITE (6,160) Z0UMB,K,0MEGA(K),CT(K),X 
rtRITE (6,150) 

FORMAT ('O', T9,'* DENOTES ORIGINAL FREQUENCY POINTS') 

FORMAT (8X, Ai , I 3,5X/2G10.9,5 a, 2R10.a,5X,2G10.9) 

IF (N.NE.3) RETURN 
WRITE (6,180) KSYM,IT£R 

format ('O', 25X, 'SYSTEM NO. ',13,', ITERATION NO. ',IA) 

DO 270 1=1 ,NM 
K=KPTS(I) 

IF ( I.EQ.KMIN + I ) GO TO 190 
IF CI.EQ.I) GO TO 220 
60 TO 2il0 
v'iRITE (6,200) 

format ( '0',25X, 'ATTEFUATED FREQUENCY INFORMATION'/) 

WRITE (6,210) 

FORMAT ('0',T2, 'NO. ',T7, 'MARGIN RAD lUS ', T26 , 'FREQUENCY ', Ta I , 'DESIR 
lED MARGIN', T57, 'MARGIN T YPE ' r T70 , 'ACTIVE '/) 

GO TO 290 
WRITE (6,230) 

WRITE (6,210) 

FORMAT ('O', 25X, 'RELATIVE STABILITY INFORMATION'/) 

XDuMb=STBM(I) 

YDUMB=RQ(I) 

IF (TYPE(I).NE.XXPO GO TO 250 
XDUMB=RA02*AS1N(XDUMB/2. ) 

YDUMB = RAD2*ASlN(Y0iJMB/2.) 

v(RITE (6,260) I,X 0 UMB, 0 MEGA (K) ,YDUMB,TYPE(I) , ACTIVE(I) 

FORMAT (' ',T2, l2,T8,GI0.9,T20,Gl0.a,T31,G10.9,TA3,G10.4,T63,Al,.T7 
12»A3) 

CONTINUE 

RETURN 


STABILITY INFORMATION'/) 


END 
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SUBROUTINE PARTAL 


SUBROUTINE PARTAL (OMEGA ,NFREO»CT#KPTS/TYPE/T»P»G»KSYM,Cl»CI r WORK I 

1 ,A,B,C,0,E,F,Z,PFX,PPY,KIN,K0UT,ZA,ZB,2C,ZD,ZE,PA,PB,PC,PD,PE,N1 , 
22,M1 ,M2,GAIN,K0NT) 

COMPLEX ANUM,PCG,G(A,B,D) ,0MEGAC0),Q,XV,CT(D) ,0IST,PGXC3) f TCAf A,D) 
1 ,P'(B,A,D) 

INTEGER A,B,C,0,E,F,Z 
INTEGER TYPEfXXT 

dimension N1(A,B), N2(A,B), M2(A,85, ZA(A,B,C), ZBCA,B,CJ 

1, ZC(A,B,C), Z0(A,B,C), ZECA,fi,C), PA(A,8,CD, PBCA,B,C)/ PC(A,8,C) 

2, PO(A,B,C), PE(A,B,C), K0NT(A,6), GA1N(A,B) 

DIMENSION KPTSCU, TYPECl), XXT(«), PFX(E,Z), PFY(E,Z) 

DATA IBLANK,XXT /HH , 1 HG , 1 HP , 1 HS, 1H A/ 

NOP=0 

DO 190 k= 1,NFREQ 
K/iHlCH=KPTS(K) 

SGN=+1 . 

IF (TYPECKKEO.XXT(il)] SGN=-l. 

IP (TYPE(Kl.EO.XXT(in GO TO 10 
IF (TYPE CK) .EO. XXT (2n GO TO 30 
IF (TYPE(K) .EO.XXT(3n 0=CMPt. X (-1 . , 0 , ) 

IF (TYPE(K).EO.XXT(i»)) 0=CMPLX(0.,0.) 

GO TO 60 
DO 20 L=0,1 

IF (AIMAG(CTCKNHICH + Ln*AIMAG(CT(KWHICH+L-l)).LE.O.) GO TO 50 
DO ao L=0,1 

IF (CCA8SCCTCKwHICH + L)I-l.)*{CA8S(CT(KWHICH+L-in-l.).LE.O.) GO TO 
1 50 

DIST=CT(KNHlCH+L)-CT(i<rtHlCH+L-l) 

XV = CONJG(CT(Kif,HICHD +1 . ) 
nistsCONJGCOISn/CARSCOISTl 
DIST=CMPLX( AIMAG(DIST) ,R£AL(DIST)) 

IF (REAL (DIST*XV).GT.0.D OISTs-OIST 
O=CTCK.hhICH)+ 5.*0IST 
DIST=CO;'iJG(-’Q + CT(KlrtHICH)) 

xv=omega(kwhichi 

khoT=0 

LNOT=0 

DO 180 1 = 1, KIN 
DO 180 J=l,KOUT 
IOP=KONT{I,J) 

CALL CRT (2,KWHICH,KSYM,T(1/1,KWHICH) ,KIN,K0UT,CT,PC1 ,1 ,KWHICH) ,PC 
IG, I, J,C1 ,CI,vvORKI,A,B,0) 

PCG=PCG*GAIN(I, J) 

NCOMO=N1 (I,J) 

IF (NCOMO.EO.O) GO TO 90 




00 80 N=lrNCOMO 

IF (N.GT.n iop=a 

aNUM=Z^-(I, J»N)+ZB(I,J,fJ)*XV 
PGX(1)=G(I-, J/KwriICH)/.ANUM 
PGXC2I=PGX(1)*XV 
DO 70 L=l,2 

70 PFX(K,KWOT+U=2.*REAL(PCG*PGX(U*OIST*^SGN) 

IF C10P.EQ.2) PFX(K,KNOT+i)=0.0 
80 K'gOT=KNOT+2 

90 NC0MD=N2CI, J) 

IF CNCOMD.EQ.O) go to 120 
DO 110 N=l,'gCOHD 
IF (N,GT.U I0P=2 

AMU‘^=ZCCI,J,N] + (ZOCI,J,NUZE(I, J,N)*XV)*XV 

PGXCnsGU, J,K/JHICH)/6NUM 

PGX(2)=PGX{n*XV 

PGXC3)=PGX(2)*XV 

DO 100 L=l,3 

100 PFXCK,KfJ0T + U=2.*REAL(PCG*PGXCL)*0IST*SGN) 

IF noP.EG.2) PFX(K,KNOT+1)=0.0 
110 KNOT=KNOT+3 

C OEflOMirjATOR PARTIALS 
120 NC0*^D=M1 (I,J) 

IF (NCOKD.EO.Ol GO TO 150 • 

DO lUO N=l,WCOMD 
IF (N.GT.n IOP=2 
Afgu’'‘-=PA(I,J,N)+P8(I, J,N)*XV 
PGXd )=GCI/J»KrtHICH)/ANUM 
PGX(2I=PGX(1)*XV 

IF nOP.EQ.2) PGXfl )=CMPLX(0.,0.) 
on 130 L=l,2 

130 PFy(k,LM0T+U=-2.*REALCPCG*PGX(L)*DIST*SGN) 
IF CIOP.EQ.2I PFY(K,LMOT + n = 0.0 

lao lnot=lmot +2 

150 NCOMD=M2CI,J) 

IF (ImCOMD.EO.OI go to 180 
no 170 N=l,NCOMO 
IF CN.GT.l) I0Ps2 

AMU’^=PCCI,J»N)t(PD(I,J,N)+PE(I,J»N)*XVjAXV 
PGX( 1 )=G(I, J,KtTlHICHI/ANUM 
PGXC2)=PGX(1)*XV 
PGX(3)=PGXC2)*XV 

IF CIOP,E0.2) PGXCn=CMpLX(0,,0.) 

DO 160 L=l,3 

160 PFY(K,LN0T+L)=-2.*REAL(PCG*PGX{L5*0IST*SGN) 
IF (IOP.EO.2I PFYCK,LMOT+1)=0,0 
170 Lf'IOT=LNOT + 3 
180 COfJTINUE 
190 CONTINUE 
RETURN 
C 

END 
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SUBROUTINE PHASEM 


subroutine PHASEM (GT0TAL»KP0INT,NM,FQM1N,F0MAX,KPTS,STBM, OMEGA) 

C SUBPROGHAV FOR CALCULATING PHASE MARGINS 

:'C 

C DEFINITIONS OF I/O VARIABLES 
C 

C GTOTAL-COMPLEX array of COMPENSATED OPEN LOOP FREQ, RESPONSE 
;C KPOINT-NO. OF POINTS' 

C OMEGA -ARRAY OF FREQS, 

C NM -COUNTER 

C KPTS -FREQUENCY NOS. WHERE MARGINS OCCUR 
C STBM -stability MARGINS OF MARGINS 

C FQMIN -LOWER FREQ. FOR MARGIN DETECTION 

C FOMAX - UPPER FREQ, FOR MARGIN DETECTION 
C 

DIMENSION GTOTALCI ) /KPTS (1) ^STBM(n» OMEGA (!)■ 
complex omega, GTOTAL 
C ■ 

P=1.0 ■ 

I DO 5 I=1,KP0INT 

S0=CA5SCGT0TAL(I) ) 

s^=sn-i .0 
IFCI.EQ. nsi=S2 

IF(AIMAG(.OMEGA(n ) .GT.FQMAX)PETURN 
IFCAIMAGCOMEGACn KLT.FQMINIGO TO 2 
IF(AbS(S2) ,LT, 1 ,0E-20)GO TO 1 
SGN=S2/ABS(S2) 

IF(SI*SGN.GT,0,0)G0 TO 2 

1 11=1-1 

IF(A6S(S2) .LT.A8S(S1))I1=I 

!vM=nm + 1 ■ 

KPTS(NM) = I1' 

IFCSl .EQ.S2)Sl=l.E-09 
FRAC=S1/(S1-S2) ■ 

STBM(NM)=CABSCP + FRAC*GTOTALa) + (l ,-FRAC)*GT0TAL(I-1)) 

2 S1=S2 

3 CONTINUE 
RETURN 
END 


SUBROUTINE POLEV 


SUBROUTINE POLE V (FW, K , X , F) 

C PROGRAM FOR EVALUATING A POLYNOMIAL AT A COMPLEX FREQUENCY 
C 

C DEFINITIONS OF 1/0 VARIABLES 
C 

C Fw -VECTOR polynomial COEFFICIENTS 

C K -ORDER OF POLYNOMIAL 

C X -COMPLEX FREQUENCY OF EVALUATION 
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complex x,sum»f 
dimension FW(l) 

SUM=FrtCK+l ) 
IF(K.EQ.O) GO TO 20 
DO 10 I=K,l,-l 
SUK=F'W(I) + SUM*X 
F = SUM 
RETURN 


SUBROUTINE SRMINS 


SUBROUTINE SRMINS (GTOTALf KPOINT , NM, P, Nr FOMIN, FQMAX/ KPTS » STBM, OME 
1A5 


c Subprogram for determining the minmuns of the open loop frequency 

RESPONSE A-ITH respect TO THE -1 POINT WHEN GIVEN POINTS ON THE 
OPEN LOOP REOUENCY RESPONSE 

C DESCRIPTION OF I/O VARIABLES 

kPOINT - NUMBER POINTS OF THE OPEN LOOP FREQ. RESPONSE GIVEN 
C ClMESA - CORRESPONDING FREQUENCIES OF CHOSEN POINTS 

C KPTS -FREQUENCY NCS, aHERE MARGINS OCCUR 
C FQMIN -minimum FRQ. CONSIDERED 

-P -POINT A.R.T. A MAX. OR MIN. IS DESIRED 

C N -DETER'-'INES whether A MAX. OR MIN. IS DETERMINED 

complex omega, GTOTAL 

DIMENSION GTOTALCn, KPTS(l), STBMd), OMEGA(l) 

ASNlsO.O 
31=0.0 

IF (N.LT.O) S1=1.0E15 
DO 50 1=1, KPOINT 

IF (AIMiG{OMEGA(D) .GT.FOMAXI RETURN 
IF (AIMAGCOMEGAd)} .LT.FQMINI GO TO 50 
S2=CABSCP+GT0TAL(I))**2 
ASN2=S2-S1 

IF (ASN2*(0 10,50,10 
10 IF (ASN1*ASN2) 20,«0,II0 

20 IF (ASN1*N) 3O,il0,iJO 

30 NM=hm+1 

11 = 1-1 
KPTS(NM)=li 

STBM(NM,)=SQRTCS1) 
no Sl=S2 

ASN1=ASN2 
50 CONTINUE 

RETURN 



c 


SUBROUTINE XCHECK 


SUBROUTINE XCHECK (K IN/ KOUT , ZA , ZB , ZC , ZO , ZE , Nl , N2, DV , ZET A , PG, NAM, K 
lE/NRTRl ,NRTR2,LPV, A,8,C/0,E,F,KKK) 

C 

INTEGER A,B,C/D,E,F 

DIMENSION ZACA/B/C), ZB(A,B,C3, ZCCA/S/C)/ ZD(A,B,C)/ ZE(A,B/C)/ 
nCA,8), N2(A,B), PGCE/F)/ DV(1)/ KJUMPC3) 

C 

DO 100 I=1/KIN 
DO 100 J=l/KOUT 
K1=N1 (I,J1 

IF (Kl.EQ.O) GO TO 50 
DO «0 K=1.,K1 

IF (NRTRI.EQ.U go to 40 
IF CZA( I , J/K) .GT. 1 .E -05) GO TO 20 
IF (ZA(I,J,K).LT,0.) GO TO 110 
IF (DVCKKK+l ) .GE.O. ) GO TO 20 
KREsl 

DO 10 L=1,NAM 

10 PG(L/KKK + n=0. 

LPV=LPV-1 

20 CONTINUE 

IF CZP(I,J,K).GT.l.E-05) GO TO 40 
IF (ZBU , J/K) .LT.O. ) GO TO 110 
IF COV(KKK + 2) .GE.O. ) GO TO 40 
KRE = 1 

DO 50 L=i,NAM 
50 PGa/KKK+2)sO, 

LPV=LPV-l 
40 KKKsKKK+2 

50 K2SN2CI/J) 

IF CK2.EO.O) GO TO 100 
DO 90 K=1,K2 

IF (NRTR2.EQ. n GO TO 90 ’ 
riNsSQRTeZCCI/J.Kl/ZE(I,J,K)) 

ZTsZDd/ J,K)/(2.*WN*ZE(I/J,K)) 

IF CZT.GE. ZETA+1 .E-03) GO TO 90 
IF CZT.LT.ZETA) GO TO 110 

P>nN=(ZE Cl , J,K)*OV(KKK + l)-ZC(I, J,K)*DV(KKK+3) )/(2.*WN*ZE(I/ J,K)**2 

AINC=-STEP/PMG 

DEL=-AINC 

VkN = SORTCZCCI, J,K)/ZECIf J/K)) 

ZET = Z0CI,J,K)/C2.*v;n*ZECI/J,K)) 

IF CZET.G£.ZETA+i,E-03) GO TO 50 
IF CZET.LT.ZETA) -GO TO 40 • 

GO TO 50 

40 DEL=0El./2. 

AINC=AINC+OEL*SGN 




A0 = ZCtI,J,K)+A-INC*DVCK2 + l) 

A1=Z0CI, J»K)+AINC*DVCK2+2) 

A2=ZECI, J,K)+AIf^C*DV(K2+5) 
ZET=A1/C2.*A2*SORT(AO/A2)) 

; IF (ZET.LT.ZETA) SGN=-1,0 

IF (ZET.GE.ZETA+0.999E-05) SGN=1.0. 

IF (ZET.LT.ZETA) GO TO AO 
IF (ZET.GE.ZETAt0.999E-05) GO TO «0 
i call SELECT (STEP,STPN£W,AINC*PMG) 

' KRE=3 

50 K2=K2+3 

,60 CONTINUE 

RETURN 
C 

•SUBROUTINE SELECT (STEP,STPNEW» STPTRY) 

IF (STEP.GE.O.) STPNEIv=AMINHSTPNEW,STPTRY) 
IF CSTRP.LT.O.) STPNEW=AMAXHSTPNEW,STPTRY) 
RETURN 
C 

ENO 


C SUBROUTINE YCHECK 


SUBROUTINE YCHECK (K I W, KOUT ,N1 , N2 , ZA , ZB , ZC , ZD, ZE , KRE , STEP, PMG , ZET 
1 ,NRTH1 ,NRTR2,DV, A,e,C,D,E,K2,STPNEW) 

INTEGER A,B,C,D,E 

DIMENSION ZA(A,P,C), ZB(A,B,C), ZC(A,B,C), 2D(A,B,C), ZE(A,B,CX,- 
11(A,B), N2(A,B), DV(1) 

|C 

DO 60 1=1, KIN 
DO 60 J=1,K0UT 
K1=N1(I,J) 

IF CK1..EQ.0) GO TO 30 
DO 20 KsIhKI 

IF (NRTRl.EQ.n GO TO 20 

IF CZACI, J,K) .GE.0.3 GO TO JO 

KRE=3 

DEL=2Aa, J,K)/DV(K2 + n 

call SELECT CSTEPiSTPNEW, (DEL•^,OOOO.Ol)*PMG) 

10 IF (ZBCI, J,K) .GE.O. ) GO TO 20 

KRE=5 

0EL=^8(I,J,K)/DV(K2•^2) 

call select CSTEP,STPNEW, (DEL-»-,00000l)*PMG) 

20 K2=K2-t2 

30 CONTINUE 

K1=N2(1,J) 

IF (KUEO.O) GO TO 60 
00 50 K=l,Kl 

IF (NRTR2.Eq,ij go to 50 
SGN=1.0 

AINC=-ST£P/PMG 

DEL=-AINC 

WN=SORT(ZC(I,J,K)/ZE(I,J,K)) 

ZET=ZDCI, J,K)/C2.*WN*ZE(I, J,K)) 




IF (ZET.GE.ZETA+I.E-03) GO TO 50 
IF (ZET.LT.2ETA) GO TO tiO 
GO TO 50 

ao 0EL=DEL/2. 

ainc=ainc+deu*sgn 

AO = ZC (I. J,K)+AINC*DV(K2*-1) 
Al=ZDCl,J,K)+AINC*DV(K2f2) 
A2=ZE(l,J,K)+AINC*DV(K2f3) 
ZET=A1/(2.*A2*SQRTCA0/A2)) 

IF (2ET.lt. ZETA) SGN=-1.0 
IF (ZET.GE.ZETA+0.999E-03) SGN=t.O 
IF (2ET.LT. ZETA) GO TO «0 
IF (ZET.GE.ZETA+0.999E-03) GO TO 40 
CALL SELECT CSTEP,STPNE^', AINC*PKG) . 
t<RE=3 

50 k2=K2+3 

60 CONTIMUE 

RETURN 
C 

SUBROUTINE SELECT ( STEP » STPNEW, STPTRY) 

IF (STEP.GE.O.) STPNEwsAMIM (STPNEWfSTPTRY) 
IF (STEP.LT.O.) STPNE1 v=AMAX1 CSTPNEW,STPTRY) 
RETURN 
C 

END- 
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